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Owen Phillips’s textbook, Flow and Reactions in Permeable Rocks, published in 
1991, became a classic in the field of geological fluid dynamics. This book is its 
long-awaited successor. In the intervening years, significant advances have been 
made in our understanding of subterranean flow, especially through the vast amount 
of research into underground storage of nuclear waste and aquifer pollution. This 
new book integrates and extends these modem ideas and techniques and applies 
them to the physics and chemistry of sub-surface flows in water-saturated, sandy, 
and rocky media. It describes essential scientific concepts and tools for hydrologists 
and public health ecologists concerned with present-day flow and transport, and also 
for geologists who interpret present-day patterns of mineralization in terms of fluid 
flow in the distant past. The book is ideal for graduate students and professionals 
in hydrology, water resources, and aqueous geochemistry. 

Owen M. Phillips is a Fellow of the Royal Society and a member of the US 
National Academy of Engineering. He has held academic posts at Cambridge Uni- 
versity and the Johns Hopkins University. He was awarded the Sverdrup Gold 
Medal of the American Meteorological Society for his contributions to oceanogra- 
phy, and a fellowship in the American Geophysical Union for his contributions to 
geological fluid dynamics. His Last Chance Energy Book, published in 1979, antic- 
ipated the first global energy crisis of the 1980s, while his recent research has been 
on sub-surface aquifer flows, the dispersal of contaminants and flow-controlled 
reactions in rocks. He has two other publications with Cambridge University 
Press - The Dynamics of the Upper Ocean (1966), which was awarded the Adams 
Prize from Cambridge University, and Flow and Reactions in Permeable Rocks 
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Preface 



This book is concerned with the dynamics of subterranean flows in the natural 
environment, with the transport and dispersal of contaminants that they may carry 
and the chemistry of the interactions between the matrix and the fluid that perco- 
lates through it in spatially random conduits or aquifer pores. It is intended for 
anyone with a quantitative interest in the world around them, and particularly for 
professionals and graduate students in hydrology, geology and environmental sci- 
ence and engineering. Many of the basic concepts originated in the late nineteenth 
or mid twentieth centuries, but they were usually developed in a very idealized and 
simplified form because few field measurements of actual sub-surface seepage or 
flow patterns had been attempted. Only recently have extensive and detailed hydro- 
logical field measurements been undertaken and their findings contain surprises 
that are re-defining the way we view this part of the natural world and begin to 
understand how it works. 

In writing this book, I have relied heavily on the guidance and advice of many 
colleagues, friends and students who listened patiently, corrected gently and pointed 
in new directions. In particular, I must thank Lawrence Hardy, John Ferry, Jim 
Wood and Gordon Wolman, all colleagues, Robert Shedlock of the US Geological 
Survey and Emory Cleaves of the Maryland Geological Survey and my Cambridge 
colleagues, Andrew Woods and Herbert Huppert who always had something new 
to show me. I am particularly grateful to my wife for her tireless reading of the 
manuscript and her suggestions for improvement. 



O. M. P. 
Chestertown 
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Introduction 



The relatively new scientific field of Geological Fluid Mechanics is concerned with 
applying the principles of fluid mechanics to the geological sciences. It is charac- 
terized by close interaction between carefully conceived laboratory measurements 
on geological flows and theoretical analyses that interpret the results in terms of 
basic physical principles. It was given this name by Herbert Huppert in Cambridge, 
one of its leading practitioners in a company that includes Andrew Woods, also in 
Cambridge, George Veronis at Yale, Stewart Turner at ANU, Canberra, and many 
other technically powerful and imaginative scientists. This present book concen- 
trates on the part of Geological Fluid Mechanics that involves the flow of passive 
and reacting fluids through porous or fractured geological media. In our planet, both 
hot and cold aqueous fluids have flowed or seeped through sand and fractured rocks 
for eons, modifying their composition by dissolution, chemical reaction and depo- 
sition. Great crystalline formations and mineral deposits were formed by nature 
during that time and modifications continue naturally. The study of these processes 
was always an interesting intellectual challenge, but one of no particular urgency. 
Yet within a couple of lifetimes the pace of change has exploded as a result of 
human activity. 

The contamination of our aquifers, and in turn the rivers and estuaries into which 
the groundwater flows, is the result of both deliberate and inadvertent injection of 
a variety of human, agricultural and industrial wastes, but our knowledge of the 
extent of these changes is meager. How long does it take for the contaminants to 
build up, where do they go and, if we remove the source, how long will it take for 
the contamination to flush out? What happens to the effluents from coal mines and 
paper mills, that are dumped into streams? Nuclear power plants help to satisfy 
our gluttonous appetite for energy without generating the primary greenhouse gas, 
carbon dioxide, but the bargain is Faustian - high-level, long-lived radioactive 
wastes are being stored on site at the nuclear power plants that generated them. 
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The long-term disposal of these wastes requires that they be removed from human 
contact for times much longer than all of human history. 

Geological fluid dynamics is concerned with how these natural systems work, 
with their patterns of flow and chemical reaction in a variety of geological media, 
whether sandy aquifers, layered sediments or mosaics of fractured rocks. The first 
geological fluid dynamicist was probably Henry Darcy (1803-58), although he 
certainly did not think of himself as such. A very accomplished hydraulic engineer, 
he worked on the Dijon water supply for a number of years, and towards the end of 
his life, he and two assistants conducted a series of hydraulic experiments on the 
flow of water through a vertical column partially filled with siliceous sand from 
the Saone River and a flow of water from the Dijon hospital water supply. The 
volume flux of water was measured in a gauging station, the pressure difference 
across the sand bed was measured by using two mercury U-tube manometers, and 
he found a very accurate linear relation between the two. His work was published 
as an appendix to his extensive report (Darcy, 1856) on the public fountains of the 
city of Dijon. 

Darcy’s study exemplifies the three essential ingredients in the scientific explo- 
ration of the nature of this world about us. The explorer needs to have (i) a detailed 
and soundly based understanding of the basic rules, the “laws of nature” under 
which the flows operate, (ii) a continuing contact with physical reality and famil- 
iarity with the results of whatever careful experiments, observations and detailed 
measurements that have been made on these flows, and (iii) the ability to put the 
two together. Darcy obviously knew a lot about hydraulics, he performed the exper- 
iment himself and he made the critical quantitative connection between flow rate 
and pressure gradient. 

In more recent times and on a larger scale we can discern these same three 
attributes that have guided the remarkable progress during the past 50 years of 
our sister science, Meteorology. The “laws of nature” that govern the motion and 
properties of the atmosphere are essentially the same as those outlined in the next 
chapter of this book, the conservation laws of thermodynamics, of Newtonian 
continuum dynamics, of chemical reactions, etc., supplemented in the atmosphere 
by the laws of radiation. The “physical reality” is the atmosphere itself, in constant 
motion and burdened by its increasing load of carbon dioxide and other greenhouse 
gases. Measurements of atmospheric pressure and temperature have been made for 
over 200 years, but the pace increased in the 1970s when the Global Atmospheric 
Research Program (GARP) stimulated a vast increase in the systematic observation, 
measurement and monitoring of the atmosphere that continues today. This in turn 
provided a strong stimulus to the development of “super-computers” that were 
needed to handle the new floods of data and the numerical models of large-scale 
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atmospheric motions being developed by Jules Chamey and others. New techniques 
were developed at the National Center for Atmospheric Research in Boulder, 
Colorado, for airborne measurements while remote sensing systems such as the 
atmospheric radar of David Atlas at NASA were becoming able to scan for clouds, 
rain and atmospheric motions over larger and larger volumes of the atmosphere. 
Today, on the evening TV weather report we see marvelous real-time, data-based 
computer simulations of local and regional weather that are vivid, ongoing and 
generally accurate. We receive warnings of rain tomorrow and impending dangers 
such as hurricanes a thousand miles away. 

In contrast, the quantitative base of data in the geological sciences is much 
more sparse. It is probably safe to assert that more atmospheric measurements 
are made every day, than have been made on geological flows in all of recorded 
history. Subterranean flow measurements are difficult and expensive to acquire by 
drilling and the data are sometimes classified for commercial reasons. The medium 
is solid, often hard, opaque, and complex in structure and composition. Seismic 
techniques are able to delineate internal structures, but most of the information that 
we have still comes from surface exposures and patterns of seepage. As a result 
of all these factors, the quantitative measurements that we do have are extremely 
valuable but still severely limited in number and scope. Particularly notable are 
the recent measurement programs on the dispersal of tagged fluids in aquifers, 
conducted mainly by the US and state Geological Surveys and their analogs in 
Europe. These have generated a leap forward in our understanding of the structure 
of the variations in permeability in sandy substrates and the spatially random flow 
field that percolates through them. 

The relative paucity of field data on geological flows presents a mis-match 
with the power and sophistication of modem digital computers. With few excep- 
tions, numerical simulations of geological flows have little measured data input, 
or quantitative comparison between the computer output and field measurements. 
Parameters can be chosen without observational or experimental basis, but simply 
to make the output “seem reasonable,” i.e. to be in accord with preconceptions. 
Though often presented as factual, and generating their own air of reality, these 
simulations are often quite misleading, and no more than digitally precise rendi- 
tions of a mostly imaginary world. There is little doubt that a more fruitful approach 
would include the development of relatively simple models with several essential 
ingredients: (i) the powerful but often neglected physical constraints such as min- 
imum dissipation, (ii) the use of measured parameters, (iii) the pertinent physical 
and chemical balances involved and (iv) the flexibility for application to a variety 
of possible structural configurations. The results must then be evaluated critically 
by comparison with whatever laboratory or field data that does exist. 
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In the next chapter of this book, the general geometrical characteristics of perme- 
able media are described, together with the basic physical and chemical balances 
that underlie the developments that follow. Two important general theorems con- 
cerning uniqueness and minimum dissipation, dating from the nineteenth century 
and often overlooked, provide useful insights into the structure of constant den- 
sity flows in complex regions with variable matrix permeability. The first part of 
Chapter 3 is a brief summary of some “classical” porous media flows, the basic 
concepts of groundwater age and the various time scales for aquifer flows. Recent 
measurements have shown that the spatially random, horizontally isotropic perme- 
ability structure of a nominally uniform sandy aquifer is associated with highly 
anisotropic dispersive characteristics of dissolved contaminants. Dissolved solutes 
in fracture-matrix media disperse rapidly in the longitudinal direction, but much 
more slowly in transverse directions. These findings can be understood best in 
terms of the minimum dissipation constraint. Unsteady flows are also of interest. 
Pressure pulses from explosions and seismic eruptions spread rapidly at acoustic 
speeds, but the residual pressure then relaxes diffusively as interstitial gas and 
liquids present flow out of fractured porous rock in seeps or geysers, at a rate that 
diminishes in time. 

Chapter 4 describes the nature of buoyancy-driven flows from convection plumes 
to freshwater wedges and gravity currents. These flows are qualitatively different 
from uniform-density flows and characteristically possess circulation in the trans- 
port velocity field. In a given geological structure, the flow patterns are no longer 
unique, which raises the possibility of instability and spontaneous evolution of one 
flow pattern into another. The archetypical thermal instability is associated with 
the name of Rayleigh (1916) and, since then, many variations of the basic theme 
have been discovered that depend on the different rates of diffusion of heat and 
dissolved salts in permeable media. These instability processes have been found 
in laboratory measurements, in contemporary natural flows and their traces left in 
ancient rocks. 

Chapter 5 synthesizes these flow patterns with the patterns of reaction, depo- 
sition and dissolution that the flow produces when the interstitial fluids and the 
matrix interact chemically. There are three dominant flow-mediated reaction sce- 
narios, reaction fronts, gradient reactions and mixing zones, each of which has 
characteristic patterns of occurrence. In many geological scenarios, the rates of 
reaction may be limited by the rate at which the flow can deliver dissolved solutes 
to the reaction site. When dissolved contaminants in a surface aquifer are absorbed 
into, or react with, the enclosing matrix, a patch of contaminant moves consider- 
ably more slowly than does the interstitial fluid, centimeters per day, perhaps. An 
extreme situation is found when the reaction involves replacement and the solutions 
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are dilute compared with the mass per unit volume of the solid reactant. The propa- 
gation speed of a reaction front may be smaller than the interstitial fluid velocity by 
many orders of magnitude, possibly being only a few centimeters per millennium. 
Application of these concepts and results can illuminate not only the formation of 
mineral deposits in paleogeologic time, but also also the accumulation, transport 
and dispersal of dissolved contaminants in present-day aquifers. 
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The basic principles 



2.1 Pores and fractures 

The geological materials with which we are concerned usually lie at one extreme 
of the range of “porous media” encountered in nature and technology. The porosity, 
the volume fraction of connected voids that allow fluid movement, may be as large 
as 0.3 or 0.4 in a well-sorted sandbank or as small as 1% in natural calcite (Pryor, 
1973). The skeletal remains of corals that abound in tropical reefs contain myriad 
interstices on scales of up to a centimeter or so and may have a similarly large 
porosity; Figure 2.1 shows a sample from Bermuda at approximately half-scale. 
This kind of structure containing fluid conduits as well as more numerous smaller 
pores is at the high-porosity extreme of those generally encountered. Compaction by 
overlying sediments, the infilling of interstices by finer grains, and the precipitation 
of cements from solution can reduce the porosity by an order of magnitude and 
reduce the permeability, as we shall see, by three orders of magnitude or more. 

Many large pores arc also apparent in dolomite from the Latemar Massif in 
northern Italy (Figure 2.2). Calcium ions from the original calcite mineral have 
been replaced by magnesium, generating dolomite. The specific volume of the 
dolomite produced in the reaction is less by 3-13% than that of the original calcite, 
so that as the reaction proceeded, the porosity increased. 

Networks of small cracks or fractures allow fluid percolation even when the 
matrix itself is relatively impermeable. Seepage from fractures can often be dis- 
cerned in roadside rock exposures. Figure 2.3 illustrates a smaller-scale network in 
a sandstone cleavage plane, made visible by stain. Stained fluid moves relatively 
rapidly through the fracture network but spreads laterally into the matrix blocks 
only slowly. Note the progression of wider, older stains passing through the whole 
sample, from which spring shorter, narrower and newer branches. 

Fault systems also provide conduits for fluid motion. Even a cursory glance at 
many field exposures often reveals layers of quartz or other minerals apparently 
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Figure 2.1. An extremely porous limestone from a Bermuda coral reef, approxi- 
mately half-scale, containing shell fragments and many interstices with scales of 
up to a centimeter or so, courtesy of Professor L. Hardie. 




Figure 2.2. Pores in dolomite from the Latemar Massif in northern Italy. The 
blocks in the scale are 1 cm long. (Photograph courtesy of Dr. E. N. Wilson.) 



deposited along fractures in a larger matrix. Figure 2.4 shows a mosaic of small 
scale fractures that have served as pathways for fluid motion until becoming filled 
by deposition from the infiltrating solution, with subsequent fractures appearing 
later. 
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Figure 2.3. A network of plane fractures provided pathways for the flow of dyed 
fluid in a sandstone cleavage plane, which then diffused into the matrix blocks, 
courtesy of Professor L. Hardie. Approximately full scale. 




Figure 2.4. In this calcite block, previous fractures have been filled with dolomite, 
while more recent fractures remain partially open. 



2.2 Geometrical characteristics 
2.2.1 Porosity 

A number of geometrical length scales are pertinent to flow through permeable 
rocks and aquifer matrices. In aqueous solutions, the intermolecular or inter-ionic 
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distances arc of order 10 -9 m, while the scale of the smallest interstices of hydro- 
logical interest is possibly 1000 times larger, i.e. 1 micron or more. As a result, 
the fluid can be regarded as a continuum whose motion through the passages 
of the medium can be described in terms of the concepts and equations of con- 
tinuum fluid mechanics. The photographs of the previous section indicate that in 
porous limestone or partly cemented sandstone, the characteristic diameter So of 
the orifices or interstices may be as large as 10 -4 -10 -3 m, with a relative few 
even larger. In sand, the distance between the individual flow paths, the diameter 
of the interstices and the grain size are all comparable. In more consolidated rocks 
and granites, however, the size of the interstices is characteristically much smaller 
than the microscale distance Z () between them. Media may also be extensively frac- 
tured on scales from 10 -2 to 10 m (see Figure 2.3) and such fracture networks are 
potentially important flow conduits. Sedimentary deposits arc frequently bedded, 
with local variations in physical properties such as the porosity and permeability 
occurring over vertical scales that are large compared with the grain size but small 
compared with the overall thickness of the bed. Finally, there arc the macroscopic 
length scales h, which specify the thickness (or smallest dimension) of the porous 
bed as a whole, and Z, its lateral extent. 

In order to relate the overall flow behavior to the average geometrical char- 
acteristics of the rocks, we must consider carefully certain statistical aspects of 
the microscopic flow through individual pores or regions of inhomogeneity, as in 
Section 2.10 below. Our primary concerns are with flow patterns and velocities, 
with the transports of heat and chemical solutes on the scale of structural variations 
or on the macroscopic length scale h of the structure itself. Immense simplifications 
are possible when the microscale Z 0 is sufficiently small compared with h that we 
can find an intermediate local scale, the matrix averaging scale / AV , that is large 
compared with the grain or matrix block size Zo yet small compared with the scale 
h of the flow patterns that we wish to resolve. Thus we require that 

<$0 < Zo Zav Zi, (2.1) 

where, as a rule of thumb, the <$C inequality signs can be interpreted to mean “is less 
by a factor of at least 10 than.” Within a given stratum, properties of the medium 
and characteristics of the flow through it, when averaged over the volume are 
expected to vary smoothly from one averaging volume to the next. When the system 
is viewed on a macroscopic scale, it can again be regarded as a continuum, with 
the “point” properties being, in fact, local averages of this kind, functions of three 
spatial coordinates (x, y, z) and possibly time t. 

When the locally averaged properties of the medium are independent of the 
position of the averaging volume, the medium is said to be homogeneous. If all 
locally averaged properties are independent of direction, the medium is described 
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as locally isotropic. A well-mixed sand bed may be, on this averaging scale, both 
homogeneous and locally isotropic. As mentioned previously, however, sediments 
are frequently deposited in such a way that the fabric preserves a record in its 
layering of the vertical direction, though there may be no differences discernible 
in the two orthogonal horizontal directions. Seasonal variations in sedimentation 
rate may produce, on a microscopic scale, a stack of horizontal laminations that, on 
a macroscopic scale, lead to different percolation characteristics along and across 
the laminations. Similarly, when elongated or plate-shaped sedimenting particles 
tend to settle horizontally, the resulting fabric will preserve a record of the vertical 
direction at the time of deposition. In general, when the medium has one preferred 
direction but its properties arc independent of rotation about that direction, it 
can be described as locally axi-symmetrical. In spite of these caveats, a uniform 
sandbank does provide the basic prototype of a classic hydrodynamic al porous 
medium, whose essential geometry, involving a three-dimensional web of minute 
intersecting fluid pathways, is found in many porous rocks and other geologic 
media at different scales and with different detailed topologies. It is convenient 
to call these “sandbank-type” media to distinguish them from the fracture-matrix 
media described later, which obey the same basic conservation laws but whose 
geometry gives them quite different flow characteristics. 

An important characteristic of a porous rock is the void fraction. The total 
void fraction r/>x is that fraction of the total averaging volume represented by the 
interstices; the solids occupy a fraction 1 — 0 T of the whole. This can be measured 
by an examination of randomly taken thin sections; since a volumetric sample can 
be considered to be a stack of plane slices, the ratio of void area to total area 
in a typical thin section is equal to the ratio of void volume to total volume, </> T . 
Similarly, along a sampling line the average ratio of total length of the line segments 
in voids to total line length is also 0 T . However, not all of the void spaces may be 
active in fluid flow through the medium. Isolated cavities or “dead end” tubes can 
contribute to (j)-y but do not provide microscopic pathways to fluid motion. A more 
significant measure for our puiposes is the connected porosity 0, in which only 
those voids that provide connections among the averaging volumes are considered. 
In general, this quantity cannot be estimated from thin-section examination without 
some additional information or assumption about the structure of the fabric, but 
it can be measured by comparison of the mass of saturated and dried samples or, 
as we shall see in the next section, by fluid observations. In this book, the term 
“porosity” refers to the connected porosity, since this is the property of interest in 
fluid motion, though in the petroleum industry it is commonly used as a synonym 
for “void fraction.” 

Clearly, 0 < 0 t < 1. At one extreme, in a material with the geometry of Swiss 
cheese, such as pumice, all the voids are isolated so that 0 — 0 while 0 T may 
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be relatively large. In poorly cemented sandstone, there may be very few isolated 
voids and (p = 4>'\ ~ 0-2 or 0.3. In consolidated rocks, f ~ 0.02-0.1, while in an 
assembly of well-sorted grains, the porosity may be as large as 0.3 or 0.4. 

A variety of textures are of interest in different contexts, and some simple 
approximate relationships are frequently useful. If, on a microscopic scale, the 
interstices consist of a complex of intersecting convoluted pathways, the porosity 
can be represented as the number of pathways per unit area found in a thin section, 
7i, times the average area of their intersections with the thin section, which is 
proportional to <$q. Thus, 

4 > ~ 77 ($0. (2.2) 



2.2.2 Double porosity in a fracture-matrix medium 

Not all hydrological flow media have the statistical geometry of a sandbank. Frac- 
tures on many scales are ubiquitous in rock strata and even casual observation of 
roadway cuttings after wet weather show them to provide pathways of significant 
fluid flow. A fracture plane (the usual terminology) should not be considered as 
an approximately uniform gap between two parallel rock surfaces, but as a more- 
or-less planar network of intersecting, ribbon-shaped fluid pathways around areas 
of close rock contact and possible accumulations of detritus. Let X represent the 
total area of the ribbon pathways in the fracture plane per unit volume of the 
medium. By again visualizing the volume as a stack of slices, it can be seen that 
X also approximates the total length of the intersections of the ribbon pathway 
network per unit area in any slice through the medium. If the mean width of the 
fluid pathways is <5p, the fracture porosity, the fraction of the total volume that they 
occupy is 

# = XS F . (2.3) 

Although <5f is the appropriate measure of the fracture aperture for specifying 
the fracture porosity, it is shown in Section 2.4 that it is far less relevant to the 
magnitude of the mean fracture flow. For a given local pressure gradient, the flow 
velocity depends on a high power of the gap width, and consequently upon the 
distribution of fracture apertures; the fluid volume flux (velocity times gap width) 
is disproportionately larger through the wider gaps. 

It is interesting to compare typical numerical values of the fracture porosity 
(p\- with values of the matrix porosity in the blocks between the fractures. In a 
moderately consolidated porous rock, </> ~ 0.2, while for the individual pathways, 
8 0 may be 10 -6 -10 -5 m, so that from (2.2), the number n of internal pathways 
intersecting unit area is very large. Fractures produced by mechanical failure of the 
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matrix may have an average gap <5 F of 10 -5 or 10 -4 m, and with a moderate density 
of pathway intersections, A ~ 3 m -1 (length per unit area), the fracture porosity is 
only about 10 -4 , which is smaller by a factor of order 2000 than that of the matrix 
blocks. Clearly, the block porosity provides the dominant reservoir for interstitial 
fluid, though as will be found later, the cracks can provide the dominant pathways 
for flow. 



2.3 The transport velocity and mass conservation 

Consider the fluid flow in a sandbank-type medium in which the mean pore velocity 
is represented by v. The porosity or active void fraction is (/>, so that over the 
fraction 0 of a unit matrix volume, the mean flow velocity is v while over the rest, 
the solids or inactive voids, it is zero. Consequently, averaged over the whole, the 
volume flux or the fluid volume flow per unit cross-sectional area or the transport 
velocity 



u = <p\. (2.4) 

The velocity u can be interpreted as the velocity with which a fluid would be moving 
if it occupied the whole space and had the same volume transport. Its direction is 
parallel to v, but since for many geological materials r/> is numerically small, u is 
substantially smaller in magnitude than v. 

The porosity of a rock sample is frequently estimated by comparing wet and 
dry weights of a sample. It could also be measured directly by taking advantage 
of equation (2.4), especially if the identification of active and inactive voids is 
difficult. If water is forced under pressure through the sample, the volume of water 
passing through per unit cross-sectional area per unit time is the transport velocity. 
If a spot of dynamically passive dye or other marker (i.e. one that does not affect 
the flow) is introduced, it will pass through the sample on average at the mean 
interstitial speed, and the ratio of the two defines </>. 

The transport velocity u is one of the primary field variables used in this book, 
since transports, not simply of fluid volume but also of heat and chemical species, 
are of primary interest. There arc two important things to remember about it. First, 
if we are concerned with the spreading of particular elements of fluid that can be 
marked with salinity, dye, or other passive contaminants, then we are necessarily 
concerned with the interstitial or pore velocity, not the transport velocity; this has 
some profound consequences, as we shall see later. Secondly, note that because 
of the particular averaging definition above, the transport velocity does not obey 
the usual rule of vector addition. If the medium is moving, being subducted, say, 
at velocity V, the interstitial fluid is indeed moving at velocity v + V, but only 
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through the fraction 0 of a transverse plane, so that the transport velocity relative 
to a fixed reference frame is </>(v + V), or u + <pV . 



2.3.1 Mass conservation 

One physical constraint on the flow patterns through porous fabrics that can be 
expressed simply in terms of u is the conservation of total fluid mass. In a fixed 
arbitrary finite region of the fabric, the mass of interstitial fluid is 

J (ppdV, 

where (j) is the porosity, p is the interstitial fluid density, and the integral is through- 
out the volume of the region. This fluid mass may change as a result of net fluid 
flow across the bounding surface, the net mass transport outward being 

J pu • dS, 

where the integral is over the surface of the region and the element of surface area 
d S = ndS is directed outward. It may also change if water is generated internally 
at the rate pS w per unit fabric volume by chemical dehydration reactions among 
the rock constituents. The net rate of change of fluid mass in the region is therefore 

— J 4>pdV — — J pu.clS + J pS w dV 

= J {-V • (pu) + P S w }dV, 

by the divergence theorem. Since the volume is taken as fixed, the initial time 
derivative can be taken inside the integral: 

J {d{cj)p)/dt + V • (pu) - pS w }dV = 0, (2.5) 

and since the volume is also arbitrary, the integrand itself must vanish. A local 
statement of fluid mass conservation is then 

^:(</>P) + V • (pu) = pS w . (2.6) 

at 

In this statement, the rate of change of interstitial water mass in the medium, 
reduced by the divergence of the mass transport velocity field, is equal to the rate 
of generation of water by rock reactions, represented by the final source term. 
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2.3.2 The incompressibility condition 



For many purposes, the complete version of (2.6) is unnecessarily general. It can 
be rewritten in the somewhat disassembled form, 

3 </> / 3 \ 

V.u = 5 w -^--U- + u-Vjln(p/po), (2.7) 

that sorts out those physical effects included in the general form (2.6) that arc usually 
less important. The first term on the right of (2.7), the source term, is significant only 
when water is being generated by chemical reaction or absorbed into the matrix. 
The next expresses the simple geometrical fact that if the porosity decreases with 
time, ( dcf> /3 1 < 0), there is a flow divergence as fluid is expelled from the decreasing 
fraction of voids. This can be important in tectonic events, but if the porosity of the 
rock changes over geological time as cementation or dissolution occurs, the time 
scales arc long, and although the cumulative changes in r/> may be large, the rates 
of change are extremely small. During and following seismic events, significant 
compaction may occur in clays or other unconsolidated sediments and f may 
change rapidly. It has sometimes been suggested that compaction is an important 
process in driving interstitial fluid flow, but when such sediments compact, the total 
volume of fluid expelled per unit volume of matrix is just the change in porosity, 
necessarily a small fraction of unity. Natural geological fluids are generally dilute 
and significant geochemical changes require movement through the matrix of 
fluid volumes that are many times the matrix volume itself. The interstitial fluid 
velocities associated with compaction arc negligible compared with those produced 
by maintained hydraulic or thermal forcing. 

Variations in space and time of the density factors in the last term of equation 
(2.7) are also usually negligible. Variations in fluid density occurring as a result of 
temperature variations may be very important in producing variations in buoyancy 
of the interstitial fluids which give rise to convective motion, but since the fractional 
change in fluid density 8p/po is very small, its influence on the flow divergence 
is minor. In geological flow fields, the individual terms in V ■ u are of order u/l, 
where u is the transport speed and / the overall flow dimension; the last term is 
smaller by order 8p/po and is therefore negligible unless boiling occurs. Another 
exceptional circumstance arises when the interstitial fluid pressure is suddenly 
altered - for example, by a natural earthquake or artificially in a pressurized well. 
The interstitial fluid is often much more compressible than the rock matrix, and 
if the pressure is suddenly reduced by a nearby fracture, the un-fractured rock 
porosity changes little but the interstitial fluid will seek to expand, forcing a time- 
dependent, local flow, a process Sibson calls “seismic pumping,” but one that may 
be more accurately described as pressure diffusion, considered in Section 3.6 of 
this book. 
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Subject to these provisos, the terms involving variations in 0 and p will be 
ignored unless they arc centrally involved in the application at hand, and in the 
absence of dehydration reactions as well, equation (2.6) simplifies to 

V • u — 0. (2.8) 



This equation expresses the statement that the volumetric divergence of the fluid 
vanishes, and is usually called the incompressibility condition. In Cartesian coor- 
dinates, with u — (u,v, w), 



du dv d w 

Vu = aI + ^ + a7 =0 ' 



(2.9) 



An integral form of this statement is frequently useful. In a fixed region of arbi- 
trary shape, the volume integral of the divergence of the velocity field is, by the 
divergence theorem, equal to the normal efflux over the bounding surface. Thus, 
from (2.8), 



J V • u dV = 



J u • dS = 0, 



(2.10) 



and the net volume flux into or out of the region vanishes. 

A complementary conceptual approach is to view the conservation of fluid mass 
in terms of individual fluid elements. As argued above, we can usually assume that 
the porosity (p is constant in time. The balance (2.6) can be written alternatively in 
terms of the mean interstitial velocity v of (2.4) as 

-j- = Tp-+v-Vp = — pV • v + 4>~ l pSw- (2.11) 

at dt 



The Lagrangian operator d/dt = 3/dl + v • V in this equation is used widely in 
considerations of fluid flow. The first paid, d/dt, represents the rate of change in time 
at a fixed spatial point, while v • V represents the rate of change that is observed at 
a point moving with velocity v up the gradient of a spatially variable property. The 
combination, written as d/dt, is then the rate of change observed in a the moving 
frame of reference when the property itself varies both in space and time, and is 
called “the derivative following the motion.” It is formally equivalent to the total 
derivative in multivariate calculus, but its physical interpretation in this context is 
useful and important. In (2.11), the rate of change of the fluid density following 
the interstitial motion is equal to minus the divergence (i.e. the convergence) of 
fluid mass plus the rate of generation of fluid mass from dehydration reactions, if 
any. If the fluid density does not change following the motion, the left-hand side of 
(2.1 1) vanishes and if there are no dehydration reactions, the last term also is zero. 
So V • v = 0 and since the porosity </> is constant, we recover (2.8). 
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Figure 2.5. In two-dimensional flow, the fluid transport between neighboring 
streamlines f and i jr + the fluid transport u Az = Ai/r is constant, so that as 

their spacing increases, the fluid velocity decreases. 



2.3.3 The stream function 



In the special case of incompressible, two-dimensional flow, the stream function is 
a quantity that provides a graphic and quantitatively accurate image of the flow, and 
has the additional advantage in analysis of reducing the number of flow variables 
by one. In two-dimensional flow, the incompressibility condition (2.8) becomes 



3 u dv 

— + — = 0 , 
dx 3 y 

which is always satisfied by a function f{x, y, t) such that 



u = 



3i fr 
dy ' 




(2.12) 



(2.13) 



The function i // is called the stream function. Clearly, if we choose a differ- 
entiable but otherwise arbitrary function xjr, then (2.13) specifies a kinematically 
possible flow field, i.e. one that satisfies the two-dimensional incompressibility 
condition (2.12). The inverse statement is also true. The statement (2.12) is the 
condition that u8x — v8y is an exact differential, 8 \j/ , say, from which (2.13) 
follows. 

The velocity component in any direction is thus represented as the gradient of 
the stream function x/r in the orthogonal direction; in particular, in the direction 
in which x[r = const., there is no transverse velocity, so that the contours t fr = 
const, represent the directions of flow at each point. Moreover, if we choose local 
coordinates with the origin at a point P, and with the x-axis chosen to lie along the 
flow direction as in Figure 2.5, then at the neighboring point Q, a distance A z from 
P, the stream function is 
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The difference in the values of the stream function at the two points is then 
Ai jr — li A", the volume flux across the line joining the two points. Streamline 
patterns contain a great deal of useful and accurate information. In a distribution 
of contours, for each of which i/r = const., the volume flux between any pair 
of these curves remains constant. Where the streamlines are close together, the 
volume transports are concentrated and the flow speed is relatively high; where they 
are widely separated, the flow speed is low. Streamlines cannot begin or end in 
the interior of a region. In re-circulating flow, the streamlines are closed; if fluid 
enters or leaves a porous region across a bounding surface, streamlines originate 
or terminate at that surface, whereas at an impermeable boundary (with no flow 
across the surface) the stream function is constant along it. Figure 3.5 gives an 
example of the use of the stream function for the representation of flow in a simple 
surface aquifer with distributed infiltration from rainfall across the water table. 

Note that the numerical value of the stream function is arbitrary to the extent 
of an additive constant, since only differences in its value at different points have 
physical significance. The value of the stream function can be assigned (usually 
zero) for one particular streamline; the values of \l/ along the other streamlines then 
give the total volume transport between that streamline and the streamline i//= 0. 
In steady two-dimensional flow, the streamlines correspond to the transport paths 
of marked fluid elements (except for diffusion), but in unsteady flow this is not 
generally so. 

In mathematical terms, use of the stream function enables us to represent the 
two velocity components u(x, y, t), v(x, y, t) in terms of the single scalar function 
i /s(x, y, t) while satisfying the incompressibility condition automatically. We then 
have one variable fewer and one equation fewer - any such simplification is always 
welcome! 

A different stream function, the Stokes stream function, can be defined for 
axially symmetrical flow, which has this same mathematical advantage but a slightly 
different interpretation in terms of the flow pattern. In cylindrical polar coordinates 
( r , 6, z), with corresponding velocity components (n, v, w), the incompressibility 
condition is 



13 1 dv 3 w 

V • u = (ru) -\ 1 = 0. 

r dr r 3 0 dz 



(2.14) 



When the flow is axially symmetrical, 3/3 9 = 0, and this reduces to 



1 3 

r dr 



( ru ) + 



dw 

dz 



= 0 , 



(2.15) 



which is satisfied by the function iff - the Stokes stream function, such that 



1 3fs 

r dz 



1 3 jr s 
r dr 



u = 



w = — 



(2.16) 
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As before, lines x/ss = const, indicate the direction of flow, so that their pattern 
gives a good visual representation of the flow field. The total volume transport 
between the axially symmetrical surfaces i// s = \//i and 1//2 is luixf — 1 // 2 ), but the 
r~ l factors in (2.16) produce a geometrical distortion in the relationship between 
fluid velocities and the gradients of V's at different radial positions. 



2.4 Darcy’s law 
2.4.1 Hydrostatics 

In a fluid at rest, the pressure increases with depth to support the weight of the 
overlying fluid; the pressure gradient is vertically downward. For a change dz in 
depth, the pressure changes by the amount pogdz, where p () is the local fluid density 
and g the gravitational acceleration. The hydrostatic pressure gradient can therefore 
be expressed as 



Vp h = (0, 0, -pog) = -pogl (2.17) 

where 1 is a unit vector vertically upward. This remains true no matter what the 
shape of the container enclosing the fluid and no matter whether the fluid as a 
whole is over- or under-pressurized. In the present context, (2.17) holds true for 
connected regions of single phase fluids at rest inside the interstices of permeable 
rocks. 

The water table is identified by the water level in an un-pumped well open to 
the atmosphere. It is defined in this book as the surface at which the interstitial 
water pressure is equal to the mean atmospheric pressure. If the groundwater is 
not moving, the water table is horizontal. It does not generally coincide with the 
upper interface of the water-saturated region, since surface tension can draw water 
upward until the capillary suction from air-water interfaces, concave on the air-side, 
supports the weight of water above the level of the water table. J. R. Philip has made 
many important contributions on water movement in soils; specific references arc 
contained in his review (Philip, 1989). In rocks, liquid water may be immobilized 
by some degree of chemical bonding with the minerals in the matrix. Because of the 
randomness in the pore geometry of most permeable rocks, the interface between 
the air and the liquid water region can be expected to be highly irregular with 
isolated pockets of water and air in the unsaturated regions above and below. Its 
geometry is self-adjusting and over short time intervals it is presumably essentially 
static, like a raft of bubbles on water. For a clean air-water interface, the surface 
tension coefficient 7c is about 73 dynes/cm (or ergs/cm 2 ) at 15 °C, though in 
field conditions, dissolved surface-active and biological materials may reduce the 
effective surface considerably. The capillary suction, the difference between the 
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air pressure and the water pressure immediately below an interface is given by 
Laplace’s formula. 



Sp=7i G; + £)- (2 - 18) 

where R\ and R 2 are the principal radii of curvature, which in this context are of 
the order of half the pore size. The liner are the pores, the greater is the capillary 
suction and the thicker is this capillary zone. In an extremely fine matrix, its 
maximum possible height above the water table may be limited by cavitation 
when the absolute pressure of the liquid drops below the vapor pressure at that 
temperature. In the laboratory under very clean and static conditions, cavitation in 
water can be avoided over a significant range of negative absolute pressures, but 
this would not be expected in natural rocks with an abundant supply of nucleation 
sites such as shaip corners, specks of debris, etc. 



2.4.2 Interstitial flow through a uniform matrix 

If the interstitial pressure gradient V p is not hydrostatic and the weight of fluid 
per unit volume pg is not necessarily uniform, the fluid will flow at a rate deter- 
mined by —V p — pgl, the driving force per unit volume. The signs express the 
facts that fluid flow is driven down the pressure gradient and that the weight of 
fluid pg acts downwards while the unit vector is conventionally taken as verti- 
cally upward. Once the interstitial fluid is in motion, internal viscous stresses are 
generated that oppose the motion, and the resulting fluid velocity in the pores is 
determined by the balance the driving and the resistive forces. In small-molecule 
fluids such as water or aqueous solutions, the viscous stresses (force per unit area) 
arc proportional to the rates of strain in the fluid, the constant of proportional- 
ity being the molecular viscosity /x, a property of the fluid. Fluids in which the 
stress, rate-of- strain relation is linear arc called Newtonian fluids. In the interstices, 
the rates of strain depend on the local geometry but are of order v/<5, where v 
represents the interstitial fluid velocity and S the characteristic pore size, so that 
the viscous stresses are proportional to /xv/<5. In an individual fluid pathway, the 
retarding viscous force acting along the walls (per unit length of path) is this 
stress times the pore circumference, which is of order 8, i.e. jiv. The driving 
force is —(V p + pg\) times the pore area of order 8 2 over which it acts. The 
balance between them can be expressed in terms of the mean interstitial velocity 
v as 



/XV oc — <5 2 (V p + pg\). 



(2.19) 
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In terms of the transport velocity u = <p\, 



u = --(Vp + pgl), 
A 



(2.20) 



where k is called the intrinsic permeability. It is proportional to cf>8 2 with a numerical 
constant of proportionality that incorporates the complicated statistical geometry 
of the connected pore spaces in the medium. This relation is known as Darcy’s law. 

An important comment should be made at this point. The linearity of the relation 
found by Darcy (1856) between the driving forces and the transport velocity is a 
consequence of both the linear viscosity relation in a Newtonian fluid (such as air, 
water or aqueous solutions) and also the neglect of inertial effects in the pore fluid 
as it moves along its convoluted pathways. If, at some point, its trajectory has a 
radius of curvature r, the fluid inertia sets up an additional pressure gradient pv 2 /r, 
where v is the pore velocity - this provides the centripetal acceleration associated 
with the curved trajectory. Darcy’s law is accurate only when these inertial pressure 
gradients are small compared with the viscous stress gradients pv/8 2 , and since in 
general r ~ 8, this requires that 



or that the combination 



p v 

T~ 



av 

«P' 



i>8 u8 

— = — = R « 1, 
v (pv 



(2.21) 



where v = p/p is called the kinematic viscosity and R is the pore Reynolds number 
which expresses the ratio of inertial to viscous effects in the flow. A corresponding 
inequality is assumed in fracture flow. This condition is usually satisfied very 
strongly in sub-surface flow except possibly for flow through coarse gravel beds 
or in vigorous geothermal systems. For such cases, Forchheimer represented the 
additional form drag by adding to the linear drag a term quadratic in the transport 
velocity and containing an empirical form drag coefficient. Other less fortunate 
embellishments of the Darcy equation include the subtraction of a term /xV 2 u from 
the pressure gradient in (2.20), the result being known as the Brinkman equation 
(see Nield, 1984) and of the addition of a cubic drag term, giving what has been 
called the cubic Forchheimer equation. Both of these, however, seem to be of more 
mathematical than hydrological interest. 

There arc, however, situations in which the simple Darcy equation is inadequate. 
If the interstitial fluid consists of two or more separate components or phases, such 
as oil and water, or air and water, the radius of curvature of the interface can be 
comparable with the pore size 8. From Laplace’s formula (2.18), surface tension 
7c can support a pressure difference of order Tq/8 across the interface without 




2.4 Darcy’s law 



21 



motion. If n represents the number of such interfaces per unit length along the 
pores, an overall pressure gradient of order n Tq /<5 can be supported without fluid 
motion. When this threshold pressure gradient is exceeded, fluids will begin to 
move, first in the pathways with largest pores and, subsequently, as the pressure 
gradient increases, in the pathways with the next largest and so on. The distribu- 
tion of pore sizes, and not simply the mean, becomes crucial in determining the 
flow characteristics, and the overall relation between pressure gradient and trans- 
port velocity is certainly nonlinear. Interesting numerical experiments on capillary 
displacement and percolation of immiscible fluids, showing fingering, trapping, 
and incomplete displacement of one fluid by the other, have been conducted by 
Chandler, Koplik, Lerman and Willemsen (1982), but this very important though 
specialized topic will not be pursued in this book. 



2.4.3 Permeability 

The permeability is an intrinsic property of the medium, just as the viscosity is an 
intrinsic property of the interstitial fluid. The permeability k has physical dimen- 
sions L 2 and is measured in units of cm 2 or m 2 and so forth. A curious unit, but one 
still widely used in hydrology, is the darcy, defined as the permeability that allows 
a transport velocity of 1 cm/s of a fluid with viscosity 1 centipoise (10 -2 e.g.s. 
units, close to that of water) under a pressure gradient of 1 atmosphere per cen- 
timeter! Since mean atmospheric pressure is ~10 6 e.g.s. units, a permeability of 1 
darcy is about 10 -8 cm 2 , or 10 -12 nr. A related quantity also widely used when 
pressures are expressed in terms of the equivalent hydraulic head is the hydraulic 
conductivity 



Pgk _ gk_ 
p. v 



( 2 . 22 ) 



where v = p/p is called the kinematic viscosity. The hydraulic conductivity 
depends on both the permeable medium and the fluid; since its physical dimensions 
are [AW 1 j (as arc those of velocity), it is expressed in units of cm/s, m/s, or m/yr, 
etc. It is a particularly convenient quantity in nearly horizontal groundwater flow 
in unconfined aquifers in which the pressure is close to hydrostatic (see below) and 
the Darcy equation reduces to the simple two-dimensional form 



u = -kV H (, 



(2.23) 



where Vh = (d/dx, d/d y) is the horizontal gradient operator and / (x, y, t) is the 
water table configuration. In a water-saturated medium, the numerical value of K 
in m/s is about 10 7 larger than the numerical value of k in m 2 . 
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The formulation of the Darcy equation (2.20) does not involve assumptions about 
any specific geometry of the fluid pathways in a three-dimensional medium, though 
the derivation above does assume that the pathway scale can be represented by the 
characteristic pore diameter 8 alone. In the proportionality statement k oc (f>8 2 , the 
numerical values of the proportionality constant are found to be generally quite 
small. Some simple geometries can be solved explicitly. If the interstices form an 
ensemble of parallel tubes of diameter 8, the Poiseuille solution of 1840 (see, for 
example, Batchelor, 1967) leads to a permeability of k = (0<5 2 /32) cos 9, where 
6 is the angle between the axis of the tubes and the pressure gradient. With an 
isotropic distribution of tube directions, the permeability is cj)8 2 / 96, even smaller 
than the previous case because many of the tubes are transverse to the pressure 
gradient. In more complex geometries, the permeability cannot be calculated so 
simply (if at all!) but it is of interest to note that while the permeability in each 
case is proportional to </><S 2 , the constant of proportionality is small, in the range 
(1-3) x 10 -2 . If, in a particular - sample, the porosity r/> = 0.1 and the pore size 8 = 
0.1 mm, we would expect k to be in the range (1-3) x 10 -11 nr, or (1-3) x 10 -7 cnr. 
Most laboratory measurements of permeability or saturated hydraulic conductivity 
involve the measurement of volume flux through a specimen under a large imposed 
pressure gradient and the use of equation (2.20), or more recently, in a centrifuge 
where the centrifugal potential gradient provides the driving force (Nimmo and 
Mello, 1991). 

Note also that the grain size per se is irrelevant; the porosity involves the size of 
the pores. The two are proportional only in the case of well-sorted, un-cemented 
grains, needles, or other particles of approximately uniform shape and size. In 
spite of this, attempts have been made to express various sets of empirical data on 
permeability in terms of porosity and grain size (see, for example, the summary 
in Lerman 1979), but the empirical formulas obtained should be used for media 
similar to those measured or, at most, only as a very general guide. 

In a useful review (concerned mainly with crystalline rocks). Brace (1980) points 
out that laboratory measurements of permeability in small samples give values for 
sandstone of from 10 -12 to 10 _16 m 2 and for limestone-dolomite, a wide range 
from 10 -14 to as small as 10 _22 m 2 , reflecting the variations in micro-structural 
characteristics. Metamorphic and granitic rocks have permeabilities from 10 -16 to 
10 -21 nr. In situ measurements give generally higher values, 10 _9 -10 -12 nr for 
chert and reef limestone, 10 _12 -10“ 18 nr for granites or crystalline rocks, probably 
reflecting flow through fractures rather than through the rock itself. Stober (1996) 
reports si mi lar values in Black Forest granites and gneisses. Fractures may be 
sealed or absent in shale. Many limestones and sandstones are already so permeable 
(~10 -15 nr) that fractures add little. In the face of these variations, it is clearly 
rash to depend heavily on an a priori assumption or estimate of average values, 
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Table 2.1. Characteristic magnitudes of permeability and porosity for some 
common rocks, orders of magnitude only. The higher values for in situ 
measurements probably result from flow through fractures, 
rather than through the matrix itself 



Material 




Permeability (m 2 ) 




Laboratory samples 


In situ measurements 


Porosity 


Sandstone 


10~ 12 -10~ 16 


o 

1 

0 

1 

o 

1 

£ 


5 x 10~ 2 


Limestone Dolomite 


10“ 14 — 10" 22 


10" 10 — 10 — 12 


(5-20) x 10~ 2 


Metamorphic/Granitic 


10“ 16 -10~ 21 


10" 12 -10“ 16 


10 -2 



Data from Brace (1980) and Stober (1996). 



but fortunately, as we shall see in subsequent chapters, flow patterns can often 
be inferred without such numerical guesses. Estimate of flow magnitudes, which 
arc critical in many practical applications involving contaminant dispersion and 
extraction, do usually involve directly the numerical values of k and their spatial 
variability. Although field calibration is sometimes possible, this often remains a 
primary uncertainty (see Table 2.1). 



2.4.4 Reduced pressure and buoyancy 

In many groundwater flow situations, the pressure remains quite close to hydrostatic 
and in the Darcy equation (2.20), the two terms on the right are almost equal in 
magnitude but opposite in sign. Their sum (which drives the flow) is much smaller 
than the magnitude of either. It is often more accurate in calculations, and is also 
conceptually revealing in flows driven by heat or salinity, to subtract out from 
(2.20) the large term V p H = — pogl representing the vertical pressure gradient in a 
fluid whose density p 0 is the average for the interstitial fluid in the domain of flow. 
The difference between the actual or total fluid pressure p\ and the hydrostatic 
pressure referred to a convenient origin for z is called the reduced pressure or 
the non-hydrostatic pressure. Since the total pressure is of little interest in flow 
considerations (though important in thermodynamics), the unadorned symbol p 
will henceforth refer to the reduced pressure p-\ — /; H . With this understanding, 
after division by p 0 , (2.20) becomes 

u = -(— V(p/p 0 ) + b\), 
v 



(2.24) 
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where the kinematic viscosity v = fi/po and the combination 




is the interstitial fluid buoyancy, a body force per unit mass that acts vertically 
upward where the local fluid density p is less than the mean p (l . Where p > po, 
the fluid has excess weight locally and the negative buoyancy is a body force 
per unit mass acting downward. Buoyancy has the same physical dimensions as 
acceleration ( LT ~ 2 ). This form (2.24) of the Darcy equation exhibits clearly on 
its right-hand side the two major driving forces for sub-surface fluid motion. The 
first term is the gradient of reduced pressure, generally associated with variations 
in the water table level above, driving fluid flow down the gradient from higher 
to lower pressure. The second expresses the buoyancy of the interstitial fluid. If, 
because of variations in temperature or salinity in the medium (greater temperature 
or lower salinity), p < po, the fluid density is less than the mean and the buoyancy 
is positive, providing a body force upward. 

Note particularly that the gradients in any horizontal direction of reduced pres- 
sure and total pressure arc the same, since by definition the hydrostatic pressure 
gradient is vertical. 

To specify the reduced pressure p = p T — p H itself, rather than its gradient, we 
must take a reference level, say z = 0, for the hydrostatic pressure and, relative 
to this, the hydrostatic pressure p H = ~Pgz, with z taken vertically upward. If 
the water table is at z — f(x, y), where the actual pressure p T is atmospheric (i.e. 
zero), then the reduced pressure at the water table is 

P = Pg£(x,y)- (2-26) 

This is equivalent to the potential introduced by Hubbert (1940), a somewhat 
redundant concept that is used widely in the groundwater literature. Its usefulness 
is limited to constant density flows in which the velocity is proportional to the 
gradient of a scalar (the potential). In flows where buoyancy is important, the idea 
of a potential is no longer useful, but the concept of reduced pressure remains 
physically pertinent. 



2.4.5 Boundary conditions 

There are two general types of condition that must be satisfied at the boundaries 
of a flow domain in order to define the flow. So-called kinematical boundary 
conditions assert in general that fluid cannot be created nor disappear - at a boundary. 
Dynamical boundary conditions assert that, except for capillarity, the fluid pressure 
is continuous across domain boundaries. A discontinuity in the pressure field 
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would imply an infinite pressure gradient and infinite transport velocity, which is 
physically unacceptable. Thermal boundary conditions specify either the boundary 
temperature or heat flux, whichever is physically more appropriate. Some specific 
examples follow. 

(1) There can be no fluid transport across a fixed impermeable boundary, so that the 
normal component of the transport velocity vanishes at the boundary. Symbolically, 
u • n = 0 there. Streamlines cannot cross such a surface; in two dimensions the surface 
must coincide with a streamline T = const. In the absence of buoyancy effects, 
Darcy’s equation shows that the pressure gradient normal to the boundary also vanishes: 
dp/dn — 0. 

(2) At a submerged boundary, between a saturated sub-strate and a lake above, the pressure 
is hydrostatic and the reduced pressure is the same as that at the water surface above. 

(3) At internal boundaries between regions 1 and 2 of different permeability, say, the 
normal component of velocity is continuous across the interface; Ui n = m • n, where 
the direction of the unit normal n is the same on each side. Consequently, from Darcy’s 
law, k\ ( dp/dn) l = ki ( dp/dn) 2 so that the normal pressure gradient on each side is 
inversely proportional to the permeability on that side. 

(4) At all points along an internal interface, the pressure is the same on each side. Conse- 
quently the pressure gradient along the interface is the same on both sides, and from 
Darcy’s law, the tangential velocity on each side is proportional to the permeability on 
that side: Ui -t/k\ — ui ■ t/C, where t is a tangential unit vector. If the permeability is 
discontinuous at an interface between two rock types, then so is the transport velocity 
along either side. Conditions (3) and (4) can be combined to produce an analogue of 
Snell’s law in optics, k sin 0 — const, where 9 is the angle between the direction of 
flow and the normal. 

However, if the boundary is “free,” i.e. an interface between two identifiable 
bodies of fluid such as fresh and saline water, or air and water at the water table, 
it may move through the matrix. The dynamical boundary condition is again that 
the pressure is continuous across the moving interface, except for capillary effects. 
The appropriate kinematical boundary condition however, is slightly more complex. 
Consider the conservation of interstitial fluid at a free boundary z — f(x, y, t) with 
the z-coordinate vertically upward, as illustrated in Figure 2.6. If fluid is flowing 
beneath a fixed interface with slope V £ , the vertical component of the fluid velocity 
is vr • V£, where Vh is the horizontal component. If, in addition, the interface is 
moving upward with velocity £ , the total vertical velocity v z of fluid at the interface 
is 



Vz = f + Vh • V£ (2.27) 

or, in terms of the transport velocities (u, w), 

w = 4>t; + u • V£. 



(2.28) 
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Figure 2.6. At a moving internal interface between saturated freshwater and 
saltwater regions, for example, the fluid elements on either side must remain in 
contact, with relative tangential motion in the interface and equal vertical motion 
on both sides. 



Equation (2.27) can be derived more formally by considering the total time 
derivative or the derivative following the motion of z = C (*, y, t), which moves 
with the interstitial fluid velocity. From the formula for the total derivative in 
multivariate calculus, 

dz 3£ 3(1 dx 3(1 dy 

dt 3 1 dx dt dy dt ’ 

and since Vz — dz/dt, and t> H = ( dx/dt , dy /dt), we recover (2.27). 

If the “free” surface is in fact the water table, there may be the additional factor 
of infiltration from rainfall at the rate W{t), measured in terms of length per unit 
time, i.e. velocity. This is generally considered as positive even though the rain 
comes from above, and represents an additional transport to the interface. Thus 
(2.28) becomes 



W(t) + w = <M+ u • Vf . (2.29) 

Some particular cases: with a horizontal surface and no internal flow, £ = W(t) /<f>\ 
in a steady state (or over a time average) in an aquifer recharge region, the vertical 
transport velocity just below the surface 



w = —Wit) + u • V£, 



(2.30) 
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both terms on the right-hand side being negative since the horizontal flow at the 
free surface is directed down-slope. 



2.5 Mechanical energy balances 
2.5.1 Flow tubes and flow resistance 

Useful and general hydrological insights can be obtained by considering the overall 
forcing, flow, resistance and overall energetics in a thin flow tube extending from 
the water table at a recharge region through a region of possibly inhomogeneous 
permeability, to its discharge at a lower elevation. The lateral boundaries of the 
flow tube are defined by streamlines of the flow. Let s be the distance along the 
axis of the flow tube from the entry point and let A(s) be its cross-sectional area. 
The permeability can vary along the length of the flow tube, so that k = k{s), in 
general. Under steady conditions (or in the mean) the incompressibility condition 
(2.9) ensures that the volume flux through the flow tube q = u{s)A{s) is constant 
along its length. The Darcy equation (2.24) can be rearranged to express the force 
balance among the driving pressure gradient, the buoyancy force and the flow 
resistance, 



—Vp/po + bl = uu/A(.s). (2.31) 



This can be integrated along a flow tube axis from any one designated point, 1, say, 
to another, 2: 



2 



~ J V(p/p 0 ) 



• ds -(- 



2 

J b\-ds = 
1 



2 




k ! u • ds. 



(2.32) 



The first line integral involving the pressure, integrates to p$ [p\ — pi\, the 
pressure difference between the two end points. If, in particular, points 1 and 2 
are at an aquifer recharge point and discharge point respectively, this reduces to 
g(£ i — £>). the hydrostatic head difference that drives the flow, as indicated in 
Figure 2.7. The second integral, which involves the internal buoyancy distribution 
that may arise from temperature or salinity variations in the flow region, can be 
written as 

B — J bcosOds, 

where 6 is the angle between the flow direction and the upward vertical. This is 
interpreted as the net buoyancy force along the flow tube axis, which is driving 
when b and cos 9 have the same sign (i.e. positive buoyancy in an upward flow or 
negative buoyancy - excess weight - in a downward flow) or retarding when the 
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A(s) 



Figure 2.7. A flow tube with cross-section A(s) and constant volume flux dis- 
charging into a lake or stream is driven by the hydraulic head of the water table 
in the infiltration region above the water level at discharge, plus any variations in 
fluid buoyancy along the tube, and is resisted by the flow tube resistance along the 
path. 

signs are mixed. In the last term, u • ds = (q / A(s)) ds, where q is the volume flux, 
which is constant along the length of the tube even though the cross-sectional area 
may vary. Equation (2.32) then becomes 



Equation (2.33) therefore expresses the two driving forces, the pressure difference 
between the two ends of the flow tube and the net internal fluid buoyancy B, 
which balance the viscous retardation of the fluid in the medium, proportional 
to the volume flux q. The flow tube resistance has the physical dimensions of 
(length x time) -1 . Four factors influence its magnitude: the fluid viscosity, the 
matrix permeability, the path length and the flow tube area. With given hydraulic and 
thermal driving forces, sections of high permeability provide small flow resistance, 
even over long path lengths. A “retarding layer” of low permeability may separate 
aquifers above and below, yet provide only moderate flow resistance if it is relatively 
thin (a small path length) and if the cross-sectional area of the flow tube is able 
expand sufficiently, as the fluid seeps across. 

A useful quantity in many applications is the flow tube conductance C, which is 
the reciprocal of the flow tube resistance: 



Po '(Pi - Pi) + B = qR, 



(2.33) 



where the flow tube resistance along the streamline is 




(2.34) 




(2.35) 



Equation (2.33) can be rewritten as 



d = C {p 0 \P\ ~ P2) + 5} , 



(2.36) 
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which expresses the volume flux as the flow tube conductance times the sum of 
the driving forces. If there are a number of flow tubes in parallel with the same 
driving forces, the flux in each is proportional to the conductance and the total flux 
is proportional to the sum of all the flow tube conductances. Low-permeability, 
low-conductance inclusions may have very little effect on the overall volume 
flux, provided the flow tubes can detour around them, as they may through a 
gap in an otherwise confining layer. The flow simply avoids regions of the lowest 
permeability. In any flow domain, the hydrologically important regions arc those of 
high conductance, whose permeability is greatest and connectivity most extensive. 



2.5.2 Energy balances 

The mechanical energy balance in the interior of the flow can be derived in a similar 
fashion. The primary energy supply is the inflow of potential energy associated 
with infiltration of the water at higher reduced pressure (or water table elevation) 
relative to the discharge, and is given by p^ x q(p\ — pi)- There may also be internal 
energy sources or sinks associated with the buoyancy distribution. The internal 
redistribution of this energy flux can be exhibited by taking the scalar product of the 
transport velocity u with Darcy’s equation for the force balance in the form (2.31). 
From elementary calculus, u • Vp = V • (pu) — pV • u and V • u = 0 because of 
incompressibility. This leads to the mechanical energy balance throughout the flow 
region, 

— p^'V • ( pu ) + b( 1 • u) = vu 2 /k. (2.37) 

The first term in this equation specifies the spatial rate at which the potential energy 
is decreasing in the flow direction, and the second expresses the rate at which energy 
is supplied or extracted from the flow by buoyancy, according to the sign of the 
buoyancy and the direction of flow. The final term, which is always positive, is the 
rate e at which energy per unit mass of interstitial fluid is dissipated by the viscous 
flow through the interstices: 

u 2 

e = v — , (2.38) 

k 

where v is the kinematic viscosity. 

The integrated version of the mechanical energy balance (2.37) over a thin flow 
tube of cross-sectional area A(s) is also of interest. Since, again, the volume flux 
q = u{s)A{s) is constant along the flow tube and with use of (2.24), the integrated 
equation becomes 

Pq 1( 1 (Pi ~ P 2 ) + q J bco$6ds = q 2 R 



(2.39) 
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in terms of the flow tube resistance R in (2.34). In this integral energy balance, the 
rate of input of potential energy that drives the flow plus the rate of energy supply 
from buoyancy is balanced by the internal rate of energy dissipation q 2 R. Note 
that the internal pressure distribution has disappeared from this overall balance 
since it provides for the flux of energy from one point to another throughout the 
interior, but is not an internal energy source or sink. The overall energy balance is 
the integral of (2.37) over the whole flow domain (i.e. over all the flow tubes): 

—Po / pu ■ ds + J bwdV — J sdV . (2.40) 

The three terms represent, in order, the net potential energy input associated with 
infiltration into the region (u • ds < 0) at higher reduced pressures and discharge 
at lower ones, the energy supplied by the upward motion of more buoyant fluid 
in the interior, or by downward motion of denser fluid, or lost by the reverse, and 
finally, the viscous dissipation throughout the entire domain. 

In these derivations, the matrix and the fluid have been regarded as incompress- 
ible. Transient compressibility effects can be important following seismic events 
as discussed in Section 3.6, but usually, when compared with the gravitational 
potential energy, the strain energy of a hydrological system (the energy of elastic 
compression of the fluid) is a trivial paid of the overall energy balance. The strain 
energy per unit mass of water is defined as 

dV dp 

Es = ~P~rr = P — - 

V p 

where dV and dp are the changes in specific volume and density produced by 
the pressure p. Changes in pressure and density are related by dp — c 2 dp, where 
c is the speed of sound, which, for the most compressible constituent, water, is 
about 1400 m/s. In a permeable formation extending to a depth h , the fluid pressure 
is approximately pgz, where z is measured down from the water table, and the 
volumetric strain in the fluid at depth 2 is dp/p = dp/ pc 2 = gz/c 2 . Consequently 
the mean strain energy of the water (occupying the volume fraction (j) of the water 
column) is 



Es ~ \<Ppgh(gh / c 2 ) 

while the mean gravitational potential energy E P = \pgh (the numerical factor 
arising because the mean pressure is half the base value). The ratio of mean strain 
energy to mean gravitational potential energy is equal to c/rgh/c 2 . Even if the 
vertical extent h of the aquifer is as large as 1 km, this ratio is less than 1%. The 
fluid kinetic energy is far smaller still. 
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2.6 Two theorems 

This section concerns the application of two theorems that were established by 
Helmholtz (see Lamb, 1932), concerning the flow of viscous fluid of constant 
density at very low Reynolds number. For a modern discussion, see Batchelor 
(1967). In principle, Helmholtz’s general results for flow in arbitrary connected 
domains do cover interstitial flow in permeable media with arbitrary distributions 
of (positive) permeability, but direct proofs for Darcy flow are very much simpler. 
Despite their antiquity, the theorems seem to have been largely overlooked in this 
area, though they do offer the potential for many conceptual insights and practical 
applications. 



2.6.1 The uniqueness theorem 

The uniqueness theorem established below asserts that there cannot be more than 
one Darcy flow solution for constant density flow with given boundary conditions - 
the solution is unique. The immediate utility of this theorem is that, if one finds 
a solution to such a problem, there are no others. The solution cannot develop 
a bifurcation or an instability. The theorem is true no matter what the internal 
permeability distribution may be. Note, however, that this uniqueness applies to 
fluids of constant density; it does not extend to situations in which buoyancy forces 
are involved. In fact, we know by examples that it is not true in these situations. One 
simple such case involves a saturated permeable region with horizontal isopycnals 
(lines of constant density). One solution is a state of rest - the velocity is everywhere 
zero. When fresh water lies over more saline water, the density decreases with 
height and the state of rest is stable - if it is disturbed slightly, it will return to its 
initial state of rest. However, when the salinity and density increase with height, 
the state of rest is unstable; tiny perturbations amplify and the system moves to a 
new time-dependent or steady-state solution. 

To prove the uniqueness for constant density Darcy flow, consider the flow 
through a medium in which the permeability k(x) is an arbitrary (but positive) 
function of position. In the absence of buoyancy forces, 

u = -(*(x)//i)Vp, (2.41) 

where p is the reduced pressure. The incompressibility condition is V • u = 0. Let 
us suppose that the theorem is false, that two different patterns of flow and pressure 
(u, p) and (u', p') are possible in a given region V, both satisfying (2.41) and 
the incompressibility condition, with assigned distributions of either pressure or 
normal component of transport velocity on the boundaries S', 

p = p or u • dS = u r • dS on the boundary surface S. 



(2.42) 
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We prove that the two solutions are, in fact, the same. Consider the following 
integral throughout the volume 

J *(x)[Vp -Vp'fdV 
= J k(x)[Vp - Wp'] ■ [Vp - V p']dV, 

= —pc J (u — u') • V(p — p’)dV , from (2.6.1), 

— —pc J V • [(u — u ')(p — p')]dV, from incompressibility, 

= —pc J ( p — p'){ u — u ) • dS. from the divergence theorem 

= 0, (2.43) 

the last step expressing the identity of boundary conditions (2.42). The first integral 
(2.43) therefore vanishes and since k(x) is everywhere positive, it follows that the 
integrand must be zero and consequently V p = V // everywhere. The two solutions 
are identical and the theorem is true. Solutions are unique. 



2.6.2 The minimum dissipation theorem 

The minimum dissipation theorem is not only conceptually important but also use- 
ful in some practical situations by providing a means of inferring, without detailed 
calculation, the general nature of flow patterns in perhaps complex geological sit- 
uations. The theorem statement is as follows. In a region occupied by a permeable 
medium in which buoyancy effects are negligible, if the transport velocity distri- 
bution over the boundary of the region is prescribed, then the actual internal flow 
has a total rate of dissipation of energy that is less than any other conceivable, 
kinematically possible flow in the same region with the same boundary conditions. 
(A kinematically possible flow is any velocity field one might imagine that satisfies 
the incompressibility condition, but not necessarily the Darcy equation. The actual 
flow, of course, satisfies both.) For example, in a near-surface groundwater flow, 
rainwater infiltrates downward to the water table at a rate that can be regarded as 
given, moves a possibly significant distance as groundwater through a matrix with 
a complex distribution of permeability and ultimately discharges into streams or 
lakes. Whatever the internal permeability structure may be, this remarkable the- 
orem asserts that the patterns of groundwater flow speed and direction are those 
which minimize the overall dissipation rate. 

The proof is as follows. Let u(x), p(x) represent the true solution with u pre- 
scribed on the boundary S. Consider a kinematically possible alternative flow that 
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one might dream up, u + u', p + p', with the same velocities at the boundary, so 
that 



u' = 0 on S. 



(2.44) 



This alternative flow also satisfies the incompressibility condition so that V • u' = 0, 
but it may not satisfy the Darcy equation. From (2.36), the total rate of energy 
dissipation in the alternative flow is 



/ 



sdV = p 
= F 



J {(u + i i') 2 /k}dV, 

J {{u 2 + u' 2 )/k}dV + 2p J 



(u • u '/k)dV. 



Since the true solution does satisfy the Darcy equation pu/k = — V p, the last term 
above becomes 



2 J (Vp) ■ u 'dV = 2^ V • ( pu')dV 

= 2f„W 



■ dS 



= 0 , 



since V • u' = 0, 

from the divergence theorem, 
from (2.44). 



Accordingly, the total dissipation rate in the alternative flow is 
m / {(u 2 + u ,2 )/k}dV, which is greater than that of the true flow, namely 

» / <“7wv . This establishes the theorem. 

It is a very important and far-reaching result. Examination of the derivation above 
confirms that it remains true if the permeability k(x) in the region is not uniform 
provided only that the velocity distribution across the boundaries is specified. 
Regions or lenses of low permeability are generally sites of proportionately low flow 
velocities, and because of the quadratic dependence of the dissipation on u, they 
contribute relatively little to the overall dissipation. The flow occurs preferentially 
in the high-permeability, low-resistance regions, attesting to the accuracy of the 
old adage that “flow follows the path of least resistance.” It is a robust theorem, 
finding a number of applications in later sections of this book. 



2.7 The thermal energy balance 

The distribution of temperature in the fabric is constrained by an equation describing 
the heat balance in each averaging volume - an expression of the first law of 
thermodynamics. Let C represent the specific heat at constant pressure; we will use 
subscripts “M” and “F” to refer to properties of the saturated matrix as a whole and 
to those of the fluid, respectively. The rate of change in time of the heat content in 
any element of unit volume is (pCm)3T /3 1, and this may come about as a result of 
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heat addition (or subtraction) brought about by convective heat transport in moving 
interstitial fluids, by molecular conduction through the matrix, and possibly by heat 
generation (or absorption) in chemical reactions. 

The convective heat flux across unit area of the fabric is given by the product 
of the heat content per unit volume of the interstitial fluid (pC) F 7\ times the mean 
interstitial fluid velocity v with which it moves, times the fraction (j) of the area 
occupied by the fluid, that is, by (pC) F u7\ The rate of accumulation of heat in the 
volume element associated with this is the net flux inwards, or minus the net flux 
outwards across the surfaces of the element of volume, 

-(pC F ) J (Tu) ■ ndS = — (pC F ) J V • (Tu)dV, 

from the divergence theorem. The convective heat input per unit volume is conse- 
quently 

— (pC) F V • (u T) = — (pC) F u • V7\ (2.45) 

in virtue of the incompressibility condition (2.8). 

Heat is also transferred by conduction through the saturated matrix. The conduc- 
tive heat flux down the gradient is proportional to the magnitude of that gradient - 
this is Fourier’s law of heat conduction - and can be expressed as — «- M V7, where 
/cm is the thermal conductivity of the saturated medium. The rate of accumulation 
of heat per unit volume is then minus the divergence of this, or 

V • (/c M Vr) - K M V 2 7\ (2.46) 

when the medium is thermally homogeneous. Finally, if chemical reactions are 
taking place, heat may be added at the rate Q per unit volume. Combining these 
expressions, we have the thermal energy balance 

d T 

(pC)u— = — (pC) F u • VT + K M V 2 r + Q. (2.47) 

at 

This can be rearranged slightly and divided throughout by (pC) f to give 

d T 

M- hu-VT = kV 2 T + Q, (2.48) 

dt 

where M = (pC)m/(pOf ? the heat source Q = Q/(pC)y, and the thermal dif- 
fusivity k , with physical dimensions (length) 2 x (time) -1 , is the matrix thermal 
conductivity divided by (pC) F . The relative density of most crustal rocks is about 
2.6, while the specific heats are between 0.19 and 0.21, so that M is generally 
0.5 ±0.1. Small glass beads are sometimes used as the porous medium in con- 
vection experiments, and for them, M ~ 0.4. The form of the advection term in 
the heat conservation equation above indicates that in a permeable, water-saturated 
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medium, the effective advection velocity for heat is u/M, somewhat larger than the 
transport velocity but smaller than the interstitial fluid velocity u/<p. This is because 
heat is diffused (or “leaked”) from the fluid pathways into the matrix solids, while 
the interstitial fluids and inert solutes remain in the interstices. 

While the advection of heat by the mean flow is always important in thermally 
driven flows, the dispersion of heat about the mean flow streamlines by the ran- 
dom excursions of fluid among the grains seems less so, except possibly in more 
permeable media such as insulating material (Kvernvold and Tyvand, 1980). The 
thermal diffusivity for saturated rock, k is of the order 10 _7 m 2 s _l (MKS units), 
about the same as that for water. Heat is conducted through the entire fabric, not just 
along the fluid pathways, which occupy only the fraction (j) of the area across which 
thermal dispersion occurs. Molecular heat conduction through the fabric dominates 
the random heat dispersion on the pore scale / about the mean streamlines when 

k (pvl, (2.49) 

a condition that is usually satisfied in geological contexts. For example, even in a 
sandy aquifer with porosity (p ~ 0-3, pore diameters of order 10 -3 m and a mean 
interstitial fluid velocity of 1 km/yr ~3x 10 -5 m/s (which for most aquifers is a 
high value) the left-hand side of (2.49) is about 10 times larger than the right. 



2.8 Dissolved species balance 

The patterns of flow through sediments are of interest to the environmental engineer 
in questions of the dispersal of reactive chemical wastes in surface aquifers and to 
the geochemist largely because of the physical and chemical changes in the fabric 
that are being produced or were produced in the distant past. Both are influenced 
by cementation in which solid material is deposited from solution at the pore 
boundaries, by dissolution when it is removed and carried away, or by various 
chemical reactions between the solid constituents of the fabric and the solutes of 
the interstitial fluid. 

Equations for the balances of chemical constituents in the interstitial fluid will 
be developed in the same way as in the preceding section. Since our interest is 
in spatial distributions of chemical reactions and reaction products, the intersti- 
tial fluid concentrations, c(x, t), are defined as mass of dissolved solute per unit 
volume of solution, with dimensions \ML~ 3 ], rather than per unit mass, which is 
dimensionless. The solute concentration balance, like the mass balance (2.11) and 
the heat balance (2.48), has the same form regardless of the units in which it is 
expressed, whether moles or mass per unit volume, provided they are consistent 
throughout. 
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The spatial distributions of fabric alteration involve processes on a wide range 
of scales, (i) On the largest scales, from that of the formation itself to the matrix 
averaging scale (Section 2.2), solutes arc transported advectively by the mean 
interstitial fluid flow v(x, t) through the pores and fractures in the medium, (ii) 
Within the matrix averaging volumes, solute is dispersed in spatially random fluid 
pathways and possibly by small-scale layering, while finally, (iii) on the scale 
of the pores themselves, molecular - diffusion to and from the pore walls allows 
sorption or interfacial chemical reaction between the solid matrix and the pore 
fluids, as described by Compton and Unwin (1990). Only the first of these is 
resolved explicitly in the solute balance equations, with dispersion within the matrix 
averaging volumes being expressed in terms of a dispersivity and the chemistry 
being specified in terms of kinetic reaction rates and solution concentrations. 

Thus, the mass of a particular solute in the interstitial fluid per unit volume 
of the matrix is (pc(x, t), where (p is the matrix porosity, and its rate of change 
is the net result of (i) the advective flux convergence —V • {</>v(x, t)c(x, t)) = 
— V • (uc), where u = cp\ is the transport velocity, and (ii) small-scale dispersive 
flux cpDVc , where D is the macroscopic dispersion coefficient (Section 2.10, 
below). The dispersive flux divergence is then cpDV 2 c. Finally, (iii) for the present, 
the sub-pore and molecular-scale processes will be expressed as an overall source 
term cpQc, where <2c represents the mass of solute added per unit volume of 
the fluid per unit time, which depends on the chemistry, the surface properties 
of the internal fluid/solid interfaces and their microscale geometry. Accordingly, 
the species balance for each dissolved constituent is represented by 

(p^- = -V-(uc) + </>Z)V 2 c + (/>£c- (2-50) 

at 

If the rate of generation of fluid in the reaction is zero or negligible compared with 
the external fluid flux, then V • u = 0 and the equation can be written in terms of 
the mean interstitial fluid velocity as 

dc dc , 

- = — +v-Vc=DV 2 c + Q c , (2.51) 

dt at 

The advection velocity for non-reacting dissolved species is the mean interstitial 
velocity v = u /</>, and (2.51) states that for fluid elements moving with this veloc- 
ity, the rate of change following the motion of solute concentration dc/dt is the 
result of both dispersion and chemical reaction. Sorbtion and chemical reaction 
with the matrix generally reduce the effective advection velocity, as described in 
Chapter 5. The difference between the effective advection velocities for heat and 
dissolved salts is brought about by the fact that advected heat changes the temper- 
ature of the whole matrix whereas non-reacting advected salts remain in the fluid 
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phase. This difference has profound consequences in the transport of dissolved 
contaminants and in thermo-haline instabilities as described in Chapters 4 and 5. 

The salinity S of the interstitial fluid is usually defined as the total dissolved 
mass per unit mass of solution; when a number of dissolved species is present 
S = p~ l Y,c, summed over the concentrations c of the separate constituents. A 
conservation equation for 5 follows from (2.51); since the equation is linear in c, 
we have 

9 S 

+ y.VS = DV 2 S+ Q s (2.52) 

at 

where the combined source term (2s = p _1 E<2c may be nonlinearly dependent 
upon the individual concentrations. 



2.8.1 Rate-limiting steps and the solute source term 

The generic expression Q c in (2.50) and (2.5 1) may represent the addition of solute 
to the interstitial fluid by one or more of a variety of physical, chemical or biological 
processes. Dissolution of the solid phase adds solute to the interstitial fluid and 
Qc > 0, while precipitation depletes it, Qc < 0. Oxidation or reduction reactions 
or organic decomposition may provide both sources and sinks. Most of the reactions 
of interest in this book occur at the solid/liquid interfaces at the surfaces of pores or 
fractures, on scales much smaller than the averaging volume that is implied in the 
continuum representation of equations like (2.52). This suggests that, among many 
other factors, Qc is proportional to the average area of reacting surface per unit 
volume of fabric, or approximately, to the volume fraction of the reacting mineral. 
Nevertheless, the Darcy-type averaging process distributes the sources throughout 
the averaging volume and Qc becomes a field valuable, determined in detail by 
interfacial surface morphology, the chemical kinetics (Lerman, 1979, Drever, 1982, 
Helgeson, Muiphy and Aagaard, 1984, and others) together with other molecular 
processes. 

Even in a process so apparently simple as dissolution, the accepted scheme, as 
described by Compton and Unwin (1990), involves a sequence of steps including 
(a) advective and diffusive transport of the reactant from its source to the reaction 
location, ( b ) adsorption of the reactant onto the surface, (c) diffusion of the reactant 
on the surface to a reactive site, ( d ) reaction, (e) diffusion of products away from 
the immediate reaction site, (f) product desoiption, and (g) diffusive and advective 
transport of the reaction products away from the reaction site. In the context of 
reactions between a solid matrix and interstitial fluids, the processes ( a ) and (g) are 
site-specific and specified explicitly by the advective and diffusive terms in equa- 
tion (2.51). The ionic concentration of interstitial fluid is usually several orders 
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of magnitude smaller than in the solid minerals with which reaction occurs, and 
although it may take decades or centuries for interstitial aqueous solutions to pass 
through a geological formation, the time required to supply mineralogically signif- 
icant amounts of reactant is even greater by a further several orders of magnitude. 
The time required to progress through the sequence of internal steps, ( b ) to (/'), 
is determined essentially by the time scale over which the slowest or rate-limiting 
step in the ladder occurs, and this is frequently much smaller than the reactant 
transport time. The reactions occur more rapidly than does the transport of reactant 
to the reaction site. In this circumstance, the reaction region is concentrated into a 
relatively narrow zone or front that propagates through the region in the direction 
of the mean interstitial flow and separates the region of unaltered rock ahead from 
the region behind in which the reaction has completed. This important process is 
described in quantitative detail in Section 5.4. 

The speed at which the reaction front propagates is proportional to the supply 
of reactant, but the internal structure of the front involves the sequence of internal 
steps, ( b) to (f), listed above. Within these, the rate-limiting step seems frequently 
to be determined by the reaction kinetics (a function of the activities involved) 
at activated sites on the surface of the solid mineral, so that the active interstitial 
surface area per unit volume of the fabric is also an important factor. In the present 
context, molecular diffusive exchange between the pore or fracture volume and the 
adjacent surface is rapid, taking only seconds or minutes in a 10 micron pore, and 
only seldom may be limiting. The reaction progress can represented in terms of 
an overall rate constant y , which lumps together these factors and the temperature 
dependence (and consequently varies widely), multiplied by a function of the 
relative saturation c/ce, where c E is the equilibrium concentration: 

Qc = yc E /(c/c E ) (2.53) 

where y has dimensions (time) -1 and the function /is dimensionless. The form 
of the function /can be determined experimentally; it may be dependent on pH or 
the concentration of other ions present. Over restricted ranges of the argument, it 
can sometimes be represented as a power-law relationship, and the power involved 
is called the order of the reaction. For example, the experiments of Plummer and 
Wigley (1976) show that calcite in water at 25 °C and 1 atmosphere partial pressure 
of COt dissolves as a second-order reaction as the pH increases from 3.9 when 
c = 0, to 5.9 as the calcite dissolves. 

There arc a few simple constraints on the form of the function/that determine the 
reaction progress. When a solid is dissolving, the rate of dissolution, proportional 
to / is expected to be a maximum when the surrounding solution is extremely 
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dilute, so that we can take 

/( 0) =1, Qc = yc E when c/c E -* 0. (2.54) 

Note that this condition in fact defines y as the actual dissolution rate in a very 
dilute solution, divided by the equilibrium fluid concentration, both of which are 
measured quantities. When the solution is saturated at a given temperature and 
pressure, a state of equilibrium exists between solid and saturated solution and the 
net source vanishes: 

/(l) = 0; Q c = 0 when c - c E . (2.55) 

When a solution is supersaturated c > c E , the function / is negative, and solute 
precipitates or crystallizes from the solution. A simple form of the distribution 
function/that satisfies these conditions is 

f(c/c E )= l-(c/c E )", n > 0, (2.56) 

though it has little other chemical justification. Near equilibrium, i.e. when c — 
c E — Ac, where Ac/c E is small (2.56) reduces to 

/(c/c E ) « n(Ac/c E ) and Q c ~ ny(Ac); 

the rate of dissolution is proportional to the index n, the reaction rate and the 
under-saturation Ac. Since the overall rate “constant” y is a lumped parameter in 
geological applications, its value can vary widely even with the same mineralogy. 
Lerman (1979) quotes values of y from laboratory measurements on silicate min- 
erals from 2 to 20 yr -1 or (0.6-6) x 10 -7 s -1 , whereas those inferred from deep 
ocean cores in siliceous sediments are several orders of magnitude smaller! 

In the species concentration balance, (2.5 1), the source or chemical reaction term 
is of general order yc since the distribution term is of order unity,/(c/c E ) ~ 1 . The 
advection term is of order vc/l where v is the mean interstitial fluid velocity and 
/, the scale of valuation of the fluid concentration in the flow direction. The ratio 

yl 

Da = Y —, (2.57) 

v 

known in chemical engineering as the Damkohler number, is the ratio of the time 
(l/v) taken for the interstitial fluid to move a distance / through the matrix, to 
the time scale for chemical reaction y~ l . In a reaction front, these two quan- 
tities are equal so that Da ~ 1 and l ~ v/y characterizes the thickness of the 
front. 

During the 1990s it became increasingly apparent that unicellular microorgan- 
isms, metal-reducing bacteria are common in groundwater and sediments, and are 
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extremely efficient catalysts in a variety of geochemical reactions. Lovley et al. 
(1989) showed that even under anoxic conditions, these organisms catalysed the 
reduction of Fe 3+ ions and the oxidation of aromatic hydrocarbons in groundwater. 
Geologic Fe 3+ oxides showed a large range of reducibility, approaching 100% in 
some materials (Zachara and Rodin, 1996), and Russell et al. (2003) claim that 
hydrolysis rate enhancements of order 10 6 can be achieved. A brief but useful 
review has been given by Newman (2003). The mechanisms of electron exchange 
arc not yet entirely clear, and the range of materials in which they occur is some- 
what uncertain. Flowever, when reactions can be catalyzed effectively in this way, 
the kinetic reaction rates may be so rapid that the reaction progress is limited by 
other steps in the sequence. The principal consequence of this catalysis may well 
be a widespread reduction in the thickness of reaction fronts with their propagation 
rate unchanged, since the front propagation speed is dependent on the supply of 
reactant, not the speed at which the reaction occurs. 

Somewhat paradoxically, one geochemical reaction scenario in which the reac- 
tions are clearly free of flow control involves the hot springs in places like Yel- 
lowstone in Wyoming, Rotorua in New Zealand and Pamukkale in Anatolia. Here, 
warm mineral-laden aqueous solutions bubble to the surface, depositing calcite in 
pools and terraces at an easily observable rate. The species balance equation in 
the spring waters is precisely identical to the balance (2.5 1) for interstitial fluids in 
the porosity limit <p —*■ 1 , i.e. where the solid matrix disappears and the transport 
velocity and the mean velocity of fluid elements become identical. Interstitial fluid 
velocities in pores or fractures may be a few m/yr while the fluid velocities u = v 
in the springs can be seen to be characteristically a few cm/s, greater by a factor 
of about 3 x 10 9 . The length scales of the reaction zones in the two cases may 
generally be comparable, so that the time scales involved in reactant transport in 
the hot springs are much smaller that in aquifer interstices and the rates are much 
larger, by factors of order 10 9 -10 1G . In spite of the generally warmer temperature 
in the hot springs, the rate-limiting step in the reaction sequence there is very likely 
to lie in the kinetics, not the flow. 



2.8.2 First-order reactions 

Near equilibrium, one might expect that the rate of addition of solute to the intersti- 
tial fluid by dissolution is linearly proportional to the local concentration difference 
from the equilibrium value: 

Qc = y(c e - c), (2.58) 

which also follows from (2.56) with n = 1. Reactions in which the rate of pro- 
duction or disappearance of a dissolved species is linearly proportional to the 
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difference between the local concentration and the equilibrium value arc described 
as first order. 

Dissolution and precipitation may occur simultaneously when polymorphic 
forms of a solid phase have different solubilities. In the case of silica, for example, 
the solubility of amoiphous silica (opal) is considerably greater than that of quartz, 
those of other polymorphs being intermediate (see, for example, Lerman 1979, 
p. 387). The silica concentration in the interstitial fluids is then described by 

— + 0 _1 u • Vc = DV 2 c + ^ y n (c„E ~ c), 



where y„ is the rate constant and c „ e the saturation concentration of the nth poly- 
morph. Alternatively, 



b <p 'u • Vc — DV 2 c + y(C — c ), 

dt 

where 




(2.59) 



(2.60) 



In a steady state with no flow, the interstitial fluid concentration stabilizes at c = C, 
which is the arithmetical average of the saturation concentrations of the individual 
polymorphs, weighted by their respective rate constants. Those for which c n e > C 
are dissolving while polymorphs for which c„e < C arc being precipitated. 



2.9 Equations of state 

The density of the interstitial fluid, usually an aqueous solution, depends on the 
temperature T, pressure p, and the concentrations c„ (mass per unit volume) of 
dissolved solids. The salinity S (dissolved mass per unit mass) is defined as the 
sum of the solute concentrations, 

5 = p-'sc„, (2.61) 

n 

where p is the solution density. The relationship among these quantities, 



p = p(p,T,S), (2.62) 

is the equation of state. It is generally nonlinear and, for some fluids, it has not 
been measured with great precision. The behavior of fresh or low-salinity water 
in the range 0— 10 °C is particularly anomalous. As the temperature decreases 
from about 10 °C, water contracts and increases in density until at 4 °C a density 
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maximum is attained. Below this, the water expands and its density decreases as 
the temperature falls towards 0 °C. Ice forms first at the water surface. 

In geological flow situations, variations in interstitial fluid density arc most sig- 
nificant in their influence on the fluid buoyancy - horizontal variations in buoyancy 
necessarily lead to fluid motion, as will be seen in the next chapter. When the fluid 
temperature and pressure arc sufficiently far removed from the critical and freezing 
points, and the temperature and salinity ranges in the field of flow arc sufficiently 
small, the equation of state can be represented with reasonable accuracy as a linear 
function of temperature and/or salinity, such as 

p = Po (l-aT + pS), (2.63) 

where p = po(p) is the reference density at the ambient pressure and 

a = -p~ l dp/dT, p = pQ l dp/dS (2.64) 

arc the coefficient of volumetric expansion and density coefficient for salinity, 
respectively. For saline water in the range 10-20 °C, a ~ 1.5 x 10 _4 (°C) _1 and 
p ~ 0.78 when the salinity is expressed as the dissolved mass of salt per unit mass 
of solution (Weast, 1972). 

However, when significant variations in interstitial density are produced by 
dissolution of a major fabric component or by ex-solution, the concentration may 
be close to saturation at the local ambient temperature: c = ce(T). The equation of 
state can then be expressed as 



p = p 0 (l - «s T), 



where 



as = - 



1 d {p{T , c e (T))} 
Po dT 



(2.65) 



can be described as the thermal expansion coefficient of the continuously saturated 
solution. Its magnitude and sign depend on the chemical nature of the solute. From 
(2.65) 



1 | dp dp 3ce 1 

7o\dT + d^Jf (’ 



and if the saturated solute (rather than any other that may be present) dominates 
the variations in density. 
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Consequently, the thermal expansion coefficient for a continuously saturated 
solution is 



a s = a 



ft dc E (T) 
Po 3 T 



(2.66) 



Representative values of a, ft and 3 ce/ 3T can be found in the International 
Critical Tables (Washburn, 1929) or the Handbook of Chemistry and Physics 
(Weast, 1972). Solutes that saturate at very low concentrations have as ~ a, the 
ordinary coefficient of expansion for water. For most solutions, dc^/dT > 0, 
though CaoStTt • 2FFO (natural gypsum) is an exception with ce decreasing from 
2.4 x 10 -2 g/cm 3 at 20 °C to about 2.2 x 10 -3 at 80 °C. Those whose saturation 
concentration is large and increases rapidly with temperature may have as < 0; as 
the temperature increases, so does the density, since the volumetric expansion is 
more than offset by the increased mass of solute in the concentrated solution. 

The fluid viscosity fi is also a function of temperature and, to a lesser extent, 
of pressure and salinity, and a detailed collation is given by Kestin, Khalifa and 
Correia (1981). For example, the viscosity of water decreases by a factor of about 
three from 0 °C (1.79 cp) to 40 °C (0.65 cp). However, the viscosity (appealing in 
Darcy’s law) always occurs in combination with the permeability, whose variability 
in natural formations is usually much greater than the variations in fluid viscosity. 



2.10 Dispersion 

The spatially random fluid motion in a porous medium with variable or random 
permeability is much more tightly constrained than is turbulent motion in a river 
or lake. In an Eulerian frame of reference, fixed with respect to the matrix, the 
velocity field is essentially steady in time, although it may vary on seasonal or 
climatic time scales. The constraint of uniqueness insists that, in the absence of 
buoyancy effects, there is only one flow solution that cannot evolve in time to 
another. The minimum dissipation theorem severely restricts the meandering of the 
flow in a randomly permeable medium. Only when it offers lower overall dissipation 
can a streamline deflect from the shortest path, in order to detour through a region 
of higher permeability or to deflect around one of lower permeability. There are no 
eddies in groundwater flow patterns. 

Nevertheless, it is intuitively evident that as marked, interstitial fluid elements 
move through an aquifer with randomly varying permeability or an extensive frac- 
ture network, they will tend to disperse about their mean pathways. The velocity 
field is spatially random but approximately steady in time. In a sufficiently large 
and generally homogeneous domain, variations in local permeability produce dis- 
tributions that tend to become Gaussian (Fickian) with the root mean square spread 
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about the mean streamline that is asymptotically proportional to ( Dt) 1 ^ 2 , where D 
is the appropriate diffusivity and t is the elapsed time. This asymptotic state can 
be attained with an effective diffusivity proportional to the product of the length 
scale of the variations in permeability and the consequent fractional variations in 
interstitial fluid velocity, provided the flow domain size is much larger than the 
characteristic scale of the permeability variations. 

The structural geology has a profound influence on dispersal. In a fracture- 
matrix medium, fluid in the fractures moves and disperses more rapidly than fluid 
in the blocks, so that a pattern of double dispersion emerges. When the dispersal 
is a consequence of a more-or-less random distribution of more permeable lenses 
that concentrate and redistribute the flow, the effective pathway intersection scale 
may be comparable with the flow domain size, and the ( Dt )*/ 2 asymptotic state 
is not attained at all. In these circumstances, the dispersion is scale dependent, 
with a root mean square particle displacement from the mean that increases in 
time at rates between t 1/2 and t. An understanding of dispersion when the scale 
of inhomogeneity is not very small compared with the flow domain size, requires 
much more detailed specification of the large-scale permeability structure than is 
needed above. Important contributions to understanding the relationships between 
dispersion and the statistics of large-scale anisotropy have been made by Dagan 
(1982, 1984), who called the phenomenon “mega-dispersion.” 



2.10.1 Kinematics of dispersion 

Consider the dispersal of individual elements of fluid or contaminant as they 
move through the spatially random permeability variations in a sandy aquifer 
or the intersecting conduits in an extensively fractured medium. Two separate 
questions are involved. The first concerns only the flow geometry or kinemat- 
ics: what particular characteristics of the spatially random interstitial velocity 
field determine this dispersal? This is considered below. The second question 
involves the specific dynamical interactions between the flow and the random dis- 
tribution of permeability as measured by Hess, Wolf and Celia (1992), and others. 
What arc the velocity variations produced by the largely horizontal mean flow 
moving through this distribution? This is discussed in Section 3.3, below. 

The basic Lagrangian analysis of fluid dispersal follows that of G. I. Taylor 
(1921) on “diffusion by continuous movements” in turbulent flow, where the loca- 
tions jc of individual fluid elements in a small averaging volume were traced in 
terms of their initial positions a and the elapsed time, t. The same approach can 
be applied to the dispersion of an injected pulse of marked or contaminated fluid 
by the spatially random but time independent velocity field in a natural geological 
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y( a, 0 = x - x 



O 



Figure 2.8. In the Lagrangian specification of fluid motion, each individual fluid 
element is tagged by its initial position a relative to a fixed origin and its motion 
expressed as the time derivative of its subsequent displacement x = x(a, t ). 

structure. As a marked fluid element moves through the medium, its trajectory is 
a function of its initial position and the elapsed time, x = x(a, t), as illustrated in 
Figure 2.8. Note that the initial instant, t = 0, the fluid element is at a, so that 
x(a, 0) = a. The interstitial velocity of the fluid element is the time derivative of its 
position: 



and this is also a function of initial position and time as we follow it along. Equation 
(2.67) can be expressed equivalently in a useful integral form: 



where t' in the integral is the elapsed time. This gives the displacement of the fluid 
element from its initial position as an integral of its velocity along its possibly 
convoluted trajectory. This and other fluid elements in the same averaging volume 
(which on the macroscopic resolution scale are at the same point) have trajectories 
that diverge in time as some move one way around the local inhomogeneities and 
others, another. From a given averaging volume, then, there is an ensemble of 
trajectories, which is equivalent to the result of marking the entire fluid in the 
volume and following the subsequent history of the cloud. 

The mean streamline through the point a is defined by the mean displacement of 
the centroid of the marked fluid from its initial position over a fixed time interval 
t, and averaged over the cloud. 




(2.67) 




(2.68) 



x(a, t) — a = v(a, t')df 





= vf, 



(2.69) 
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where the averaging process is denoted by the over-bar. The dispersion of a fluid 
element about the mean streamline passing the locality a is found by subtracting 
(2.69) from (2.68) 



x(a, t) — x(a, t) — f {v(a, t') — \]dt' = f v'(a, t')dt f , 

Jo Jo 

where v' is the valuation in velocity of the fluid element about its mean velocity as 
it moves along its trajectory. More concisely, 

y(a, t) — f v'(a, t')dt' , (2.70) 

Jo 

giving the displacement y of the fluid element relative to the mean position of the 
ensemble (the cloud) in terms of its velocity v' relative to the mean. The simplest 
measure of its dispersion is the second moment of the distribution, or the mean 
square value of y, and this can be expressed usefully in terms of the geometrical 
and flow characteristics. 

In aquifer flows with random spatial variations in permeability, the direction 
of the mean streamline through a point is in the direction of the local pressure 
gradient, and the dispersal characteristics in the longitudinal, transverse and ver- 
tical direction may well be different (Hess et al., 1992). In the longitudinal, 1- 
direction, say, the displacement of the fluid element relative to the centroid of the 
cloud is 



Jt(a, t) = 




v\ (a, t')dt' , 



or equivalently, 

— yi(a, t) = r/j(a, t), with yi(a, 0) = 0. (2.71) 

When the two left- and right-hand sides of these equations arc multiplied together, 
there results 

y, lT = IT,* = l ”'><». 

Upon averaging over many fluid elements from the same initial resolution volume, 
we have 



(2.72) 
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Figure 2.9. The Lagrangian autocorrelation function /l(t) decreases from unity 
when r = 0 and its integral is assumed to converge to the integral time scale 7f , 
which is defined as the width of the rectangle of equal area. 



where r is the time interval between the two velocity variations (relative to the 
mean) in the covariance of the integrand (Figure 2.9) and aj = y \ is the variance 
or mean square displacement in the ^-direction about the centroid of the marked 
fluid. 

The velocity fluctuations of a marked fluid element are statistically steady in 
time, so that this covariance depends on the time interval r, rather than on t and t' 
separately. When r = 0, the values of v[ in the covariance arc taken simultaneously, 
and it reduces simply to the mean square of the interstitial velocity component in the 
longitudinal direction, v[ 2 . As the time interval increases, the interstitial velocity 
variation at the end of the interval begins to differ in a random way from what 
it was at the beginning as each element moves through the pores or along the 
interstices. The velocity variations relative to the mean become less correlated - 
they lose memory - and the covariance function decreases. Ultimately, when the 
interval r - t is large, the velocity variation at the later time is assumed to have 
lost all dependence on what it happened to be initially. If this is so, the covariance 
function drops to zero as r — ► oo, in the manner shown qualitatively in Figure 2.9 
and the integral is assumed to converge. 

The covariance can be expressed equivalently as 



v'i(t)v[(t + r) = v\ 2 f L (r) (2.73) 

where / L (r) is the dimensionless longitudinal Lagrangian correlation function 
(following the fluid elements) which has the value of unity when r = 0 (the 
fluctuation correlates perfectly with itself) and is assumed to decrease to zero (no 
correlation) as the time interval becomes very large. The integral time scale T L , the 
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area under the curve in Figure 2.9, is a useful measure of the time interval over 
which the correlation persists. 

The rate of increase of the longitudinal variance of the cloud of fluid elements 
in the x-direction is therefore given from (2.71) and (2.72) as 

= 2v 'i 2 [ Mr)dr 
at Jo 

/“OO 

— > 2i/j 2 / / l (t)c?t when the elapsed time r » If, 

Jo 

= 2v[ 2 T l , (2.74) 



where 



71 = 



-/ 



fh(r )dr 



is the integral time scale of the longitudinal velocity fluctuations following the 
motion. Thus 



yr 



(2v[ 2 T L y = 2D x t, 



say, 



(2.75) 



and D x — (2v' x 2 X T ) is the longitudinal macroscopic dispersion coefficient, with 
physical dimensions [L 2 T -1 ]. Asymptotically, the variance in length of the marked 
cloud increases linearly in time, and its average root mean square spread increases 
as t */ 2 , factors common to most dispersive phenomena. The variances in lateral and 
vertical spread of marked fluid can be expressed in the same form, and involve the 
mean square lateral and vertical interstitial fluid velocities and the corresponding 
integral time scales. 

Under natural field conditions the locally averaged permeability may vary con- 
tinuously on the scale of the flow domain itself, producing corresponding velocity 
correlations that, albeit small, persist over considerable spatial intervals. Inte- 
grals based on measured data may not converge and the classical description 
above, though conceptually important, is less useful. An alternative approach using 
structure functions or “variograms” that arc insensitive to large-scale variations in 
medium properties is described in Section 3.3. 



2.10.2 Dispersion in a steady plume 

In many applications the spatial distribution of contaminants from a discharge may 
be of greater interest than their precise time of arrival, and the results of the analysis 
above can be re -phrased to address this. Note that the Lagrangian analysis above is 
asymptotically exact within its defining parameters - a homogeneous, statistically 
steady, random velocity field, however it is generated. It contains no dynamics and 
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makes no use of the defining property of interstitial fluid flow that the geometry 
of the interstices does not change in time. A parallel Eulerian version argument is 
only approximate, but does explicitly recognize this fact. 

Consider a statistically steady plume of marked fluid arising from a local injec- 
tion point in a statistically homogeneous aquifer. The total time derivative in (2.72) 
reduces to the advective rate of change v • V, and in the time interval dr, fluid 
elements have moved an average distance dr\ = V\dr, where v\ is the mean inter- 
stitial fluid velocity in the 1 -direction or mean flow direction, so that in terms of 
the spatial covariance of the local interstitial velocity variations on scales of / AV or 
greater, 



Rij( r) = v' it (a)v'j(a + r) 



(2.76) 



we have 



v\ 



d 

dx I 








vj(a)nj(a + r^dri/v = 2(iq) 1 



Xi 

J R\i{r\)dr\, 

o 



Asymptotically, 



(iVi) -> 2 ( v[ 2 /v 2 ) j Mri)dn = 2 (^ 2 /U 2 ) i_x , (2.77) 

o 



where f\ \(r\) is the dimensionless correlation function associated with the covari- 
ance (2.76) and A.n_i is the integral length scale of the velocity variations in the 
1-direction along the direction of flow as in Figure 3.5 in the next chapter. The mean 
square longitudinal dispersion of fluid elements therefore asymptotically increases 
linearly with travel distance 

(yf) -* 2 (^2/U 2 ) (A.!,.!)*! = 2 ( u\/U 2 ) (Au.Ox!, (2.78) 

where u\ represents the local variation in velocity about the mean U and x\ is the 
distance from the source. Similarly, the lateral dispersion is given by 

(yf) -2 {ui/U 2 ) (k 22 _ i ) Xl , (2.79) 

where k 22 _j is the integral length scale in the 1 -direction of the velocity variations 
in the lateral 2-direction. The vertical dispersion about the mean streamline is 
expressed similarly, with 3 replacing 2 in (2.79). Since the variations in mean 
interstitial fluid velocity are produced by spatial variations in permeability, they 
arc also proportional to the local mean velocity, and the first factor on the right of 
these expressions is purely numerical, generally smaller than unity, independent of 
the mean flow speed and dependent only on the local permeability structure. 
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The expressions 

oi d - 1 = (uj/U^Xn-i old- 2 = (u\/ U 2 ^j X22-1 old- 3 = 33— 1 

(2.80) 

have the physical dimensions of [length], and are called the longitudinal, transverse 
and vertical dispersivities. They are the analogues of the macroscopic dispersion 
coefficients such as D in equation (2.75) that specify the spreading in time of a 
tagged fluid patch. Since the velocity ratios in (2.80) depend upon the permeability 
structure but arc independent of the mean flow speed, the dispersivities are proper- 
ties of the medium, not of the flow. Thus, from (2.78) and (2.79), the asymptotic root 
mean square spread of marked fluid in longitudinal and lateral directions increases 
with distance x\ from the source as {a D _\X\) l l 2 and (a £) _ 2 Vi) 1 / 2 etc. 
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Patterns of flow 



3.1 Flow in uniform permeable media 

Typical patterns of flow in a saturated aquifer are qualitatively different from those 
in a stream or lake. In the latter, flows are governed by a balance between gravity 
and fluid inertia, they may be vortical with closed circulation paths or eddies that 
range in scale from centimeters or less in a turbulent patch to the whole length of the 
lake. In an aquifer, the flow is governed by a balance between the gravitationally 
induced pressure gradients (or the internal buoyancy) on the one hand, and the 
flow resistance of the medium on the other. Fluid makes its way from recharge 
areas to discharge. We show below that circulating flows in a horizontal plane are 
impossible and that closed vertical circulation can be driven only by buoyancy 
variations, either positive or negative. 

These flows are governed by the Darcy force balance (2.24): 

k ( pa — p\ 

u = -(— ' V(p/po) + b\), where the buoyancy b — g I ) , (3.1) 

^ \ Po J 

p is the reduced pressure, and 1 is a unit vector vertically upward. This is, in 
component form with axes (x, y, z.) and velocity components (u, v, w) with z and 
w being in the vertical direction, 



k d(p/p 0 ) 
v dx 



k d(p/po ) 
v 3 y 



k 



w = — 



v 



d(p/Po) 

dz 




(3.2) 
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where the permeability k and kinematic viscosity v are here assumed to be constant. 
The Cartesian form of the incompressibility condition V • u = 0 is 



du dv dw 

— + — + — = 0 . 

d.x dy dz 



(3.3) 



Note the structure of these equations. The buoyancy, associated with variations 
in temperature and/or salinity, is a field variable determined by the conservation 
equations (2.48) and (2.51), together with an appropriate equation of state. The 
buoyancy appears only in the vertical force balance. A local region of positive 
buoyancy drives the interstitial fluid vertically upward, yet the motion is generally 
fully three-dimensional. Coupling with the horizontal motion occurs through the 
incompressibility condition (3.3). As the buoyant region moves upward, adjacent 
fluid moves inward near the bottom to replace the fluid driven vertically and outward 
near the top. 



3.1.1 Flow constraints 

In the absence of buoyancy variations, the range of possible flow solutions in a 
uniform medium is very tightly constrained. The Darcy equation reduces to its 
classical form 




P 



Note that both the permeability k and viscosity /x = p 0 v are positive and the 
minus sign indicates that the transport velocity is always down the gradient of 
reduced pressure. An immediate consequence of this is that: 

(i) in constant density Darcy flow, no streamline can form a closed loop. 

The proof of this statement by reductio ad absurdum is simple and quite general. 
Let us assume that it is not true, that there is indeed a closed streamline, and we find 
that this produces a contradiction. Start at some point on the streamline and move 
in the direction of flow; the reduced pressure would continually decrease, and as 
we complete the circuit and return to the same point, it would be less than it was at 
the beginning. But this is not possible since the pressure is a single-valued function, 
so that the assumption of a closed streamline must be false. Consequently, there 
cannot be any closed streamlines in a constant density Darcy flow, either in two or 
three dimensions, regardless of the non-uniformity or distribution of permeability. 

A dynamical characteristic of permeable-medium flow pattern, second only in 
importance to the transport velocity is its curl, the rotation vector, 



£2 — (£2 V , £2 V , £2-.) = V x u = curl(u). 



(3.5) 
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The rotation vector has the same appearance as the vorticity in ordinary fluid 
mechanics, but because it appears in the context of Darcy flow, it does not have 
the important dynamical properties possessed by the vorticity. It is, nevertheless, 
a useful property of flows driven by buoyancy and those in media of variable 
permeability, considered later in this book. If the permeability is uniform, the curl 
of the Darcy equation (3.1) is 

k 

£2 — Vxu=-{Vx (M)} . 
v 

Since 1 is a unit vector vertically upward, the vertical component of Q vanishes, so 
that 

(ii) for Darcy flow with variable buoyancy in a uniform medium, the hori- 
zontal flow field is irrotational. 

A useful mathematical theorem due to Stokes provides an association between 
the rotation vector and transport flow circulation, as illustrated in Figure 3.1. This 
theorem (see Lighthill, 1966, pp. 55-57, or any good book on vector calculus) 
states that for any vector field u(x), the line integral around an arbitrary closed 
circuit (not a streamline circuit, because none exist) of the component of u(x) in 
the direction of the circumference, is equal to the surface integral of the normal 
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Figure 3.2. In constant density Darcy flow, the line integral around any closed 
loop of u • ds/k vanishes; in particular, if the permeability is uniform, the rotation 
vector is everywhere zero. 



component of the rotation vector £2 = V x u = curl u over any surface capping 
the circuit 






u • d s 




xu ■ dA — 



J Q-dA, 



(3.6) 



where ds is a differential element of the loop, a vector whose magnitude is the 
length of the element and whose direction is always that of the tangent to the loop 
in a consistent sense and dA is an element of area of the cap. With a small circular 
contour of radius r, the left-hand side of (3.6) is equal to the mean azimuthal speed 
u times the circumference 2: rr, and the right-hand side is equal to the mean rotation 
normal to the circuit times the area nr 2 . Thus nr 2 Q = 2nru and so Q. — 2 u/r. 
The angular velocity associated with the distribution of transport velocity is u/r, 
or, from the last equation, The rotation can thus be visualized as twice the 
local angular velocity associated with the transport velocity field. 

Further dynamical constraints can be established by considering the flow at 
points along a closed loop in a constant-density flow region, as illustrated in 
Figure 3.2. This loop cannot be entirely along a streamline in view of the proof 
above, but must cut across streamlines over at least some part of the loop. Consider 
(3.4) in the form 
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and integrate this around the closed loop. Notice from Figure 3.2 that u • ds = 
u(ds) cos 9 is positive around the top of the loop since 9 < n /2, but negative 
around the bottom where the angle 9 is obtuse, so that contributions to the integral 
around the whole loop include both positive and negative parts, which tend to 
cancel out. In fact, they cancel out exactly, even if the permeability varies along 
the path. The line integral around the complete loop is 

u ■ ds f 

—j— = - (h Vp ■ ds = pi - P 2 , 

which vanishes because the circuit is closed, the finishing point 2 being the same 
as the starting point 1, and the pressure is single-valued. Thus, 




(iii) even though the medium permeability k varies spatially (with the texture 
or composition of the matrix), in the absence of buoyancy variations 
the integral around any closed loop of the tangential velocity divided 
by the permeability vanishes: 






u • ds 
k(x) 



= 0 . 



(3.7) 



This is already a useful expression. Note that it is consistent with the flow boundary 
condition at an internal interface described in Section 2.4, i.e. that the ratio of the 
tangential velocity components on each side is equal to the ratio of the permeabil- 
ities, as the reader may verify by taking a thin loop that runs along one side of the 
boundary a distance ds, then across and back on the other side and across to the 
starting point. 

From Stokes' theorem, equation (3.7) is equivalent to 



V x 




whence, by writing out the components of the curl, 



H = Vxu = V 




X u, 



(3.8) 



(3.9) 



where Q is the rotation vector of the transport velocity field and ko is any convenient 
reference value of the permeability field. Note that if a different reference value is 
chosen, the logarithm is different by an additive constant, which is immaterial since 
the permeability appears only in the gradient of the logarithmic term. The rotation 
vector theorem (3.9) demonstrates that variations in log-permeability generate 
rotation in an otherwise uniform stream and provides a general but direct connection 
between the velocity and permeability distributions without specific reference to 
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the pressure distribution. It asserts, for example, that in a largely horizontal flow, 
streamlines cutting across a log-permeability gradient induce rotation about the 
vertical axis, attracting streamlines to regions of high permeability and repelling 
them from regions of low permeability. It is used in Section 3.3 concerning fluid 
flow dispersion in an extensive aquifer with spatially random permeability. 

In a homogeneous region of constant density and permeability, (3.7) becomes 

<j) u • ds = 0, (3.10) 

so that, from (3.6) or (3.9), all components of the rotation vector vanish, £2 = 0. 
Thus, 

(iv) constant-density Darcy flow in a homogeneous, uniform medium is 
everywhere irrotational, with the circulation around any closed loop 
vanishing, and the rotation vector everywhere zero. 

This is a stronger constraint than (ii) above. It can also be seen directly from (3.5), 
with b = 0. These statements (i)-(iv) limit substantially the range of dynamically 
possible flow patterns in uniform permeable media. 



3.1.2 Laplace’s equation 

Let us return to the incompressibility and Darcy equations to show that a similar 
direct connection can be found between the permeability and pressure fields. The 
divergence of the Darcy equation (3.4) gives 

V u = — /x _1 V • (jfeVp) - -p,-\Vk • Wp + kW 2 p) = 0, 

because the fluid is regarded as incompressible and divergence-free. On rearrange- 
ment, this becomes 



V 2 p + k~ l Vk ■ Vp = 0 
or 

V 2 p + Vln(jfc/*b) • Vp = 0, (3.11) 

where k (] is a reference value for the permeability distribution, frequently taken 
as its geometric mean. Note that here, as well as in equation (3.9), the lack of 
homogeneity in the medium is expressed by the log-permeability, rather than 
the permeability itself. This equation specifies the distribution of pressure (and 
thence, transport velocity) in flow through an aquifer with a prescribed internal 
log-permeability distribution and appropriate pressure boundary conditions. 




